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Abstract 

We provide a semi-analytic study of the small scale aspects of the power spectra of warm dark 
matter (WDM) candidates that decoupled while relativistic with arbitrary distribution functions. 
These are characterized by two widely different scales k eq ~ 0.01 (Mpc) -1 and kf s = V3k eq /2 (V^)a 
with (V^)2 <C 1 the velocity dispersion at matter radiation equality. Density perturbations evolve 
through three stages: radiation domination when the particle is relativistic and non-relativistic 
and matter domination. An early ISW effect during the first stage leads to an enhancement 
of density perturbations and a plateau in the transfer function for k < kf s . An effective fluid 
description emerges at small scales which includes the effects of free streaming in initial conditions 
and inhomogeneities. The transfer function features WDM-acoustic oscillations at scales k > 2 kf s . 
We study the power spectra for two models of sterile neutrinos with m ~ keV produced non- 
resonantly, at the QCD and EW scales respectively. The latter case yields acoustic oscillations 
on mass scales ~ 1O 8 M . Our results reveal a quasi- degeneracy between the mass, distribution 
function and decoupling temperature suggesting caveats on the constraints on the mass of a sterile 
neutrino from current WDM N-body simulations and Lyman-a forest data. A simple analytic 
interpolation of the power spectra between large and small scales and its numerical implementation 
is given. 
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I. INTRODUCTION 



In the concordance ACDM standard cosmological model dark matter (DM) is composed 
of primordial particles which are cold and collisionless[]J. In this cold dark matter (CDM) 
scenario particles feature negligible small velocity dispersion leading to a power spectrum 
that favors small scales. Structure formation proceeds in a hierarchical "bottom up" ap- 
proach: small scales become non-linear and collapse first and their merger and accretion 
leads to structure on larger scales, dense clumps that survive the merger process form satel- 
lite galaxies. 

Large scale simulations seemingly yield an over-prediction of satellite galaxies [2[ by almost 
an order of magnitude larger than the number of satellites that have been observed in Milky- 
Way sized galaxies j^-I^J. Simulations within the ACDM paradigm also yield a density profile 
in virialized (DM) halos that increases monotonically towards the centerj^UlQ] and features 
a cusp, such as the Navarro- Frenk- White (NFW) profile 0| or more general central density 
profiles p(r) ~ r" 13 with 1 < /3 < 1.50,13, 1Q|. These density profiles accurately describe 
clusters of galaxies but there is an accumulating body of observational evidence 11- 17|, 19, 20[ 
that suggest that the central regions of (DM)-dominated dwarf spheroidal satellite (dSphs) 
galaxies feature smooth cores instead of cusps as predicted by (CDM). This difference is 
known as the core-vs-cusp problem[T3]. Salucci et.al. 18] reported that the mass distribution 
of spiral disk galaxies can be best fit by a cored Burkert-type profile. 

In ref.[io| a "galaxy size" problem has been reported, where large scale simulations at 
z = 3 yield galaxies that are too small, this problem has been argued to be related to that 
of the missing dwarf galaxies. 

Thus there seems to be emerging evidence that the ACDM paradigm for structure for- 
mation may have problems at small scales |2l|. 

Warm dark matter (WDM) particles were invoked 22- 24]] as possible solutions to the 
discrepancies both in the over abundance of satellite galaxies and as a mechanism to smooth 
out the cusped density profiles predicted by (CDM) simulations into the cored profiles that 
fit the observations in (dShps). (WDM) particles feature a range of velocity dispersion 
in between the (CDM) and hot dark matter leading to free streaming scales that smooth 
out small scale features and could be consistent with core radii of the (dSphs). If the free 
streaming scale of these particles is smaller than the scale of galaxy clusters, their large scale 
structure properties are indistinguishable from (CDM) but may affect the small scale power 
spectrum [251] so as to provide an explanation of the smoother inner profiles of (dSphs), fewer 
satellites and the size of galaxies at z = 3(20|. 

Furthermore recent numerical results hint to more evidence of possible small scale dis- 
crepancies with the ACDM scenario: another over-abundance problem, the "emptiness of 
voids" [26| and the spectrum of "mini- voids" [27[ both may be explained by a WDM candi- 
date. Constraints from the luminosity function of Milky Way satellites [28J] suggest a lower 
limit of ~ IkeV for a WDM particle, a result consistent with Lyman-a[29l-l3ll|. galaxy 
power spectrum 32J and lensing observations 33| . More recently, results from the Millenium- 
II simulation 3J suggest that the ACDM scenario overpredicts the abundance of massive 
> 10 10 M Q haloes, which is corrected with a WDM candidate of m ~ 1 keV. This body of 
emerging evidence in favor of WDM as possible solutions to these potential small scale prob- 
lems of the ACDM scenario warrants deeper understanding of their small scale clustering 
properties. 

A model independent analysis suggests that dark matter particles with a mass in the keV 
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range is a suitable (WDM) candidate [351. l36f . and sterile neutrinos with masses in the ~ keV 



range are compelling (WDM) candidates [37H40I . I42l448l|. These neutrinos can decay into an 
active-like neutrino and an X-ray phot on |49l| . and recent astrophysical evidence in favor of 
a 5 keV line has been presented in ref.[50| (see also(5l|). The analysis in ref.(52| suggests 
upper mass limits for a sterile neutrino in the range ~ 6 — 10 keV. Possible direct detection 
signals of such candidates have been recently assessed in ref . [53J] . 

A property of a dark matter candidate relevant for structure formation is its distribution 



function after decoupling[54j-|57j. It depends on the production mechanism and the (quan- 
tum) kinetics of its evolution from production to decoupling. There are different production 
mechanisms of sterile neutrinos |42| . l44j . |57| , leading in general to non-thermal distri- 

bution functions. There is some tension between the X-ray [I^] and Lyman-q forest (29 -3~D 
data if sterile neutrinos are produced via the Dodelson-Widrow (DW) [37| non-resonant mix- 
ing mechanism, leading to the suggestion (52j that these may not be the dominant (DM) 
component. Constraints from the Lyman-a forest spectra are particularly important be- 
cause of its sensitivity to the suppression of the power spectrum by free-streaming in the 



linear regime [29H31J. The most recent constraints from the Lyman-a forest [30l. 131 



im- 



prove upon previous ones, but rely on the(DW)[37| model for the distribution function of 
sterile neutrinos, leaving open the possibility of evading these tight constraints with non- 
equilibrium distribution functions from other production mechanisms, such as those studied 
in refs. [H3, H . 

The gravitational clustering properties of collisionless (DM) in the linear regime are de- 
scribed by the power spectrum of gravitational perturbations. Free streaming of collisionless 
(DM) leads to a suppression of the transfer function on length scales smaller than the free 



streaming scale via Landau damping [25l. |59|, |60j. This scale is determined by the decoupling 



temperature, the particle's mass and the distribution function at decoupling [61]. 

Goals: The most accurate manner to obtain the transfer function for DM perturbations 
is to use the publicly available computer codes for cosmic microwave background (CMB) 
anisotropies 0-0, with modifications that would allow to include the different distribution 
functions of the WDM particles. These codes include baryons, radiation, neutrinos and 
DM and yield very accurate numerical results. The drawbacks in using these codes for 
WDM particles are that they do not readily yield to an understanding of what aspects of 
a distribution function influence the small scale behavior, and must be modified for the 
individual WDM candidates because their distribution functions are "hard-wired" in the 
codes. 

The goals of this article are twofold: i) to provide a semi-analytic understanding of the 
main physical processes that determine the transfer function of WDM candidates at small 
scales that entered the horizon well before matter-radiation equality for arbitrary distribution 
functions, ii) to provide a relatively simple formulation of the power spectrum that allows 
a straightforward numerical implementation, valid for arbitrary distribution functions. In 
order to achieve these goals we must necessarily invoke several approximations: a) we neglect 
the contribution from baryons, b) we also neglect anisotropic stresses resulting from the free 
streaming of ultrarelativistic standard model active neutrinos. These approximations entail 
that the results of the transfer functions will be trustworthy up to 10 — 15% accuracy. 
However, the main purpose of this work is not to obtain the WDM transfer function to a 
few percent accuracy, but to provide a semi-analytic "tool", to study the main features of 
the transfer function at small scales for a particular WDM candidate given its distribution 
function determined by the microscopic process of production and decoupling. If the transfer 
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function features important small scale properties that could potentially lead to substantial 
changes in structure formation, this would warrant more accurate study with the CMB codes 
and eventual inclusion into N-body simulations. 

In this article we study the transfer function for WDM density and gravitational perturba- 
tions by solving the collisionless Boltzmann equation in a radiation and matter dominated 
cosmology including the perturbations from the radiation fluid for arbitrary distribution 
function of the WDM particle, thus the results (within the acknowledged possible errors) 
are valid for z > 2. 

Strategy: 

WDM particles with a mass in the ~ keV range typically decouple from the primordial 
plasma when they are relativistic. For example sterile neutrinos produced non-resonantly 
via the Dodelson-Widrow (DW) (37| mechanism or by the decay of scalar or vector bosons 
(BD)[44, 45, 13, 58] decouple at the QCD or Electroweak (EW) scale respectively. Therefore 



these species decouple when they are still relativistic in the radiation dominated era and 
become non-relativistic when T w m ~ keV when the size of the comoving horizon rj < Mpc. 

Therefore we anticipate that there are three stages of evolution for density perturbations: 
I) when the particle is still relativistic, this is a radiation dominated (RD) stage, II) when 
the particle is non-relativistic but still during the (RD) era, III) when matter perturbations 
dominate the gravitational potential (the particle is non-relativistic in this era). When the 
WDM particles are relativistic, their contribution to the total radiation component is neg- 
ligible because their effective number of degrees of freedom is 1 (see below). Therefore 
during stages I) and II) the gravitational potential is completely determined by the radiation 
fluid. Our strategy is to solve the Boltzmann equation for WDM density perturbations in 
the three stages. In stages I) and II) the gravitational potential is completely determined 
by the radiation fluid and the Boltzmann equation is solved by considering the gravita- 
tional potential as a background determined by the Einstein-Boltzmann equation for the 
radiation fluid. In stage III) when matter perturbations dominate, the 00-Einstein equation 
for small scale perturbations is equivalent to Poisson's equation. The initial conditions for 
the Boltzmann equation are given deep in the (RD) era when the relevant modes are well 
outside the horizon. In this work we consider adiabatic initial conditions determined by the 
primordial perturbations seeded during inflation. The main strategy is to use the solution 
of the integration of the Boltzmann equation in a previous stage as the initial condition for 
the next stage. During stage I) suppression by free streaming is independent of the distri- 
bution function and the free streaming scale grows with the comoving horizon. However 
we find that modes that enter the horizon when the particle is relativistic with wavelengths 
up to the sound horizon are amplified via an early integrated Sachs- Wolfe effect (ISW) as 
a consequence of the time dependence of the gravitational potential produced by acoustic 
oscillations of the radiation fluid. The evolution of WDM density perturbations at the end 
of this stage determine the initial conditions for stage II) when the particle becomes non- 
relativistic but still the gravitational potential is determined by the perturbations in the 
radiation fluid. During this stage the free streaming scale depends only logarithmically on 
the comoving horizon. Whereas CDM perturbations grow logarithmically during this stage 
(Meszaros effect), WDM perturbations are suppressed by a free streaming function that 
depends on the distribution function of the decoupled WDM particle. In stage III) when 
WDM perturbations dominate the gravitational potential, density perturbations obey a self- 
consistent Boltzmann-Poisson integral equation which we analyze in a systematic expansion 
valid for small scales. 
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The main results are: 



There are two relevant scales: k eq ~ 0.01 (Mpc) 1 which is the wavevector of modes 
that enter the Hubble radius at matter-radiation equality, and the free streaming scale 

kf 



'Cfs 



where (V^)^ is the mean square root velocity dispersion of the WDM particle at 
matter-radiation equality. For a WDM candidate with m ~ keV produced non- 
resonantly and decoupling either at the electroweak or QCD scale kf s > 10 3 k eq . The 
free streaming length scale l/kf s is proportional to the distance traveled by a non- 
relativistic particle with average velocity (V^L)? from matter-radiation equality until 
today, and it also determines the size of the (comoving) horizon (conformal time) when 
the WDM particle transitions from relativistic to non-relativistic: 

v/3 

Vnr 



V2k fs ' 

This means that perturbations with k > kf s enter the horizon when the WDM particle 
is still relativistic and undergo suppression by relativistic free streaming between the 
time of horizon entry until t)nr- 

During the (RD) era acoustic oscillations in the radiation fluid determine the gravita- 
tional potential <fi. The time dependence of <fi induces an early ISW that results in an 
enhancement of the amplitude of WDM density perturbations for wavelengths larger 
than the sound horizon of the radiation fluid at r] NR , namely t)nr/ \/3, but those with 
krjNR/ v3 3> 1 are suppressed by relativistic free streaming. 

In stage III), we turn the Boltzmann-Poisson equation into a self-consistent differential 
integral equation that admits a systematic Fredholm series solution. Its leading term 
is the Born approximation and lends itself to a simple and straightforward numerical 
analysis for arbitrary distribution functions. This approximation is equivalent to a fluid 
description but with an inhomogeneity and initial conditions completely determined 
by the past history during stages I) and II) . The resulting fluid equation is a WDM 



generalization of Meszaros equation |65H67j . The solutions describe WDM acoustic 



oscillations, the suppression by free streaming is manifest in the inhomogeneity and 
initial conditions. 

In the Born approximation we obtain a semi-analytic expression for the transfer func- 
tion and compare it to the CDM case. We also provide an expression for the power 
spectra that interpolates between large and small scales and give a concise summary 
for its numerical evaluation for arbitrary distribution functions, mass and decoupling 
temperature. 

We study the transfer functions and power spectra for two different scenarios of ster- 



ile neutrinos produced non-resonantly: via the (DW) mechanism[37j and via boson 



decay [5 71. |58fl. The transfer functions are very different even for the same mass. The 
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ISW amplification of density perturbation enhances the transfer function for wave- 
lengths larger than the free streaming length and is more pronounced for the colder 
species (BD). WDM acoustic oscillations are manifest in both transfer functions for 
^ ~ 2&/ s - We analyze the main physical aspects of these oscillations and suggest 
that their amplification by non-linear gravitational collapse may lead to dumpiness 
on mass scales ~ 10 8 M for the colder species . We compare the results from the 
semi-analytic approach for (DW) sterile neutrinos with the transfer function obtained 



by numerical integration of Boltzmann codes in refs. [31|, |40|, |41| • The results of the 
Born approximation agree to < 5% with the numerical fit to the transfer function 
provided in ref.[3l[ in the region where the fit is valid. 

An important corollary of the study is a quasi- degeneracy: not only the value of the 
mass but also the detailed form of the distribution function along with the decoupling 
temperature (in the form of the number of relativistic degrees of freedom at decoupling) 
determine the transfer function. Two particles of the same mass but very different 
distribution functions and decoupling temperatures, may feature very different power 
spectra. Conversely, two WDM particles of different masses and different distribution 
functions may feature similar power spectra on a wide range of scales. 

This is studied in detail here through a comparison between sterile neutrinos produced 
by the two mechanisms (DW,BD). This result suggests a caveat in the constraints on 
the mass of the (WDM) particle from current (WDM) simulations and the Lyman-a 
forest data. 



This work differs from that in ref . [68| that analyzes (standard model active) neutrinos as 
(WDM) but in an Einstein-Desitter cosmology, in two main aspects: i) our study includes 
the (RD) era and the transition to matter domination (MD) including the time dependence 
of the gravitational potential, which is a source of an early ISW effect during stage I) and 
the history during stages I) and II), and ii) we study non-thermal distribution functions. 
The inclusion of stages I) and II) during the (RD) era also distinguishes this work from that 



in ref.[56L 57 . 



II. PRELIMINARIES 



We consider a radiation and matter dominated cosmology: 



[a + a, 



eg J 



(HI) 



where the dot stands for derivative with respect to conformal time (rj), the scale factor is 
normalized to ao = 1 today, and 



1 <eq 



1 

3229 



Introducing 



(II.2) 
(II.3) 



'-eq 
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it follows that 



leading to 



da 
drj 



tt2q 
"0 m 



<eq 



[1 + a 



= 288.46 



vTT 



(Mpc) 



(II.4) 
(II.5) 



where we have used Q m h 2 = 0.134 6^|. At matter-radiation equality we define 

9.8 x KT 3 



keq — H e q d e q 



s/2 



tt2q 
"0 m 



<cq 



Mpc 



(II.6) 



corresponding to the comoving wavevector that enters the Hubble radius at matter-radiation 
equality. Furthermore from (1II.5j) we find the comoving size of the horizon at matter- 
radiation equality 

2^-1) g ^z^ 120Mpc (II7) 



from which we obtain 



V 

a = — 



Heq&eq 



1 + ^ 



Heqdeq 



Veq 



2{V2- 



V2 
k 



eq 



During radiation domination 



and in this regime 



<C 1 



(II.8) 



(II.9) 



r] ~ a 144.23 (Mpc) . (11.10) 
During matter-radiation domination, a comoving wavevector k enters the (comoving) 



Hubble radius when k = Ha corresponding to a value of the scale factor 



1 + *(&)' 



(11.11) 



We are interested in small scale properties for perturbations with comoving wavelenghts 
100 pc < A = 2n/k < 10 Mpc corresponding to k 3> k eq . For these modes, which have 
entered the horizon during the radiation dominated era, it follows that 



v/2^«l 
k 



(11.12) 



A weakly interacting massive particle (WIMP) of mass m ~ 100 GeV that undergoes 
chemical freeze-out at T c h ~ m/20 and thermal decoupling at Td ~ 10 MeV when a^ ~ 10 - ', 
and rj d ~ 10 pc, i.e, deep in the (RD) era, is non-relativistic at decoupling. Scales < r] d where 
inside the horizon when the DM particle was still coupled to the cosmological plasma and 
acoustic oscillations of the photon fluid are imprinted on the transfer function at these very 
small scales (70]. However, larger scales were outside the horizon and their perturbations 
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are frozen, they enter the horizon after decoupling and their evolution is described by the 
collisionless Boltzmann equation. 

On the other hand, sterile neutrinos with mass m ~ keV decoupled thermally at much 
higher temperature (~ 150 MeV for (DW) j37|, ~ 100 GeV for production via scalar or vector 



boson decay [5 71. |58|). and become non-relativistic at T ~ m ~ keV, namely for a ~ 10 



-3 



In terms of conformal time, m ~ keV sterile neutrinos become non-relativistic at 

77jra~0.2Mpc( — ), (11.13) 

so that for r\ <C t]nr this (DM) candidate is relativistic and non-relativistic for r\ > r/^R. 
Therefore, for DM candidates that decoupled for temperatures Td > 10 MeV all modes of 
cosmological relevance for (comoving) scales A > 50 pc may be studied in the linear regime 
via the collisionless Boltzmann- Vlasov equation. A firmer estimate will be provided in 
section ( IITTA"j) . 

For WDM particles with m ~ keV we see from eqns. (111.121 III.6P that comoving scales 
A > 0.2 Mpc entered the horizon when the DM particle is non-relativistic, whereas smaller 
scales entered during the radiation dominated stage when the WDM particle is relativistic. 
Therefore comoving scales smaller than that of cluster of galaxies became sub-horizon during 
(RD) when the WDM particle is still relativistic. This is important because free streaming 
changes from the relativistic to the non-relativistic case: during the relativistic stage the free 
streaming length is of the order of the horizon, but much smaller during the non-relativistic 
stage (see below). 

Hence as anticipated above, there are three distinct stages of evolution of density pertur- 
bations for WDM particles with m ~ keV and scales smaller than 0.2 — 1 Mpc: 

• I) (RD), relativistic r\ < 7]nr, 

• II) (RD), non-relativistic r\ eq > r\ > t]nr, 

• III) matter domination (MD), non-relativistic for 7] > i] eq . 



III. EVOLUTION OF PERTURBATIONS: THE BOLTZMANN EQUATION 



We follow the notation of Ma and Bertschinger 7_l| (see also 72M76| ) . and consider only 
scalar perturbations in the conformal Newtonian gauge (longitudinal gauge) with a per- 
turbed metric 
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i.) ■ 



9ij = a 2 (v) 1 - 20(f,r]) 
The perturbed distribution function is given by 

f(p, X, 7]) = f (p) + Fi{p, x, V) 



(III.l) 
(III.2) 

(III.3) 



where fo(p) is the unperturbed distribution function, which after decoupling obeys the colli- 
sionless Boltzmann equation in absence of perturbations and p, x are comoving momentum 



and coordinates respectively. As discussed in ref. [55M57] the unperturbed distribution func- 
tion is of the form 

where 

y = ij!- (ni.5) 

where p is the comoving momentum and T ^ is the decoupling temperature today, 

T , d =(-Y Tcmb, (IIL6) 

with g d being the effective number of relativistic degrees of freedom at decoupling, T C mb = 
2.35 x 10~ 4 eV is the temperature of the (CMB) today, and Xi are dimensionless couplings 
or ratios of mass scales. 

Although our study will be carried out for arbitrary / , we will analyze in detail two can- 
didates for WDM: sterile neutrinos produced by the Dodelson-Widrow (DW) (non- resonant) 
mechanism for which 

AW = i?TT (IIL7) 

where (3 ~ 1CT 2 37], and sterile neutrinos produced near the electroweak scale by the decay 



of a scalar with a mass of the order of the EW scale or vector bosons (BD), which are 
abundant at temperatures near the EW scale jH3, 58|, 



M P ) = X g -H2pl ; miM = f^t^- (III.8) 



and A ~ 10~ 2 57], |58|. We will compare the results for the WDM distributions with that for 
weakly interacting massive particles (WIMPs) which freeze-out with a Maxwell-Boltzmann 
(MB) distribution, 

f,{p)=Me-i ; x = ^ (III.9) 

J-d 

where m ~ 100 GeV, Td ~ 10 MeV is the thermal decoupling temperature, and M is 
determined at chemical freeze-out [781] . 

An important observation for WDM candidates is that during the radiation dominated 
era when these are relativistic, their contribution to the energy density is 

for sterile neutrinos produced by the Dodelson-Widrow (DW) or scalar decay (BD) mecha- 
nisms. Namely these WDM candidates contribute to the radiation component with an effec- 
tive number of degrees of freedom proportional to 0, A, ~ 10~ 2 and can be safely neglected in 
their contribution to the radiation component. The same argument justifies neglecting the 
anisotropic stress (quadrupole moment) arising from the free streaming of these particles 
when they are relativistic. 
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Introducing spatial Fourier transforms in terms of comoving momenta k (we keep the 
same notation for the spatial Fourier transform of perturbations), the linearized Boltzmann 
equation for perturbations is given 

by 0473 



^ftp;,,) + ^ftpi,,) + (^) [p - jfc^fM) f (M)] =0 (III.ll) 

where p = k ■ p, dots stand for derivative with respect to conformal time r\ and 



e{p,rj) = \/p 2 + m 2 a 2 {rf) 



(111.12) 



is the conformal energy of the particle of mass m. Durin g (R D) and (MD), the 00 component 
of Einstein's equation in conformal Newtonian gauge is[72| 



(j){k, ?7) + ( ^ <P(k, rj) + ^{k, rj) 



3 k 2 q 



4 k 2 a 2 



« i — ) +(!e. 

p / m \ p J r 



where 



n = - = aH 

a 



[1 + 5] 
V2a 



is the inverse comoving Hubble radius, and 



1 



5pj{k,Tl) = -7 



d 3 p 



t(p,V) F i,j{k,P,r]) 5 3 =r,m. 



a-* J (2tt) 3 

In what follows we neglect stress anisotropies leading to 

<j>(k,r)) = if)(k,rj) , 



(111.13) 



(III. 14) 



(111.15) 



(111.16) 



thereby neglecting the quadrupole moment from relativistic standard model (active) neutri- 
nos. We also neglect the baryonic component in the matter contribution, a compromise that 
allows us to pursue a semi-analytic understanding of the (DM) transfer function at small 
scales. The remaining Einstein's equations are not necessary for the discussion that follows. 
In absence of stress anisotropy, Einstein's equation (1111.131) can be written in another useful 
form, 



2k 2 a 2 



3 k 2 q 



[1 + a) (a « 



8p 



P / m 



5_p 

p Jr 



where 



d_ 

da 



(111.17) 
(111.18) 



The formal solution of the Boltzmann equation ( 1III.11I) is 



F 1 (k,p;r]) = F^k.p-rn) e 



-ik /j,l(p,ri,rii) 



-P 



dfojp) 
dp 



d<p(k, t) 
dr 



k p 



V(p,r 
(111.19) 



<l>(k,T) 
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where 

FT) 

1( P ,V,V')= / V(p,r)dr ; V(p,r) = (111.20) 

is the comoving free streaming distance that a particle travels between rf and rj with physical 
velocity V(p, r) = p/e(p, r). 

The solution ( 1III.19P with dIII.20[) is the starting point of our analysis. The density and 
gravitational perturbations produced by a WDM particle with m ~ keV that decouples from 
the plasma when it is still relativistic are obtained by evolving the solution 19[) through 
the three stages : I) when the DM particle is still relativistic during (RD), II) when the 
particle becomes non-relativistic for a > 1CT 3 but still during (RD), III) during (MD) a > 1 
(the DM particle is non-relativistic). 

During the first two stages the perturbation in the gravitational potential <p in ( 1HL19I) is 
completely determined by the radiation component to which the WDM candidate contributes 
negligibly as discussed above. The difference between stages I) and II) is manifest in the 
free streaming distance l(p, rj, rj'). During stage III) the gravitational potential is determined 
by the DM density perturbations self-consistently through Poisson's equation (this is the 
advantage of the conformal Newtonian gauge). Our strategy is to determine initial conditions 
deep in the radiation era when the cosmologically relevant modes are still superhorizon, and 
to evolve the solution f ]III.19j) through each of these stages, using the distribution function 
at the end of each stage as the initial condition for the next stage, thereby propagating the 
initial condition determined deep in the radiation era to matter-radiation equality. 



A. Free streaming distance: 



The free streaming distance l(p, rj, rf) can be obtained analytically with (jII8p . the general 
result can be expressed in terms of elliptic functions, however it is unyielding and not very 
illuminating. It simplifies considerably in two relevant cases: for radiation domination when 
V "C f]eq which includes the era when the DM candidate becomes non-relativistic, and in the 
non-relativistic regime for rj ^> t]nr which includes the matter dominated era. 

Radiation domination (RD): 

Since fo(p) is a function of y = p/T ^ d it is convenient to write p = yT 0td in V{p). In the 
radiation dominated era rj <C rj eq during which a(r]) ~ 77/77* we find 



kl(p,v,v') = -y ln 



z + 



lz' + 



V 2 a 2 _|_ z 2 



+ z 



z = krj 



(HI.21) 



where we have introduced 



fK kT ,d 

a = 2V2 — - — ! — ~ 

171 K, e q(X e q 



(HI.22) 



Since for the WDM distributions under consideration y f(y) is strongly peaked at y 
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y 2 where 



r 105 C(7) 
12 C(5) 



— = J™y 4 fo(y)dy 



15 



C(5) 
C(3) 



i 8.505 ; for (BD) 

12.939 ; for (DW or thermal fermion) 



(111.23) 



3 X 3 rri i 



for (MB) 



it follows that for z,z' <C \jy 2 oi the ultrarelativistic approximation v(p,rj) ~ 1 is valid 1 , 
and in this regime 

i(p,v,v') = (v-v'), (III.24) 

which is the comoving free streaming distance traveled by an ultrarelativistic particle be- 
tween 7] and rf . In the opposite limit when the particle is non-relativistic but still in the 

radiation dominated era z, z' \ y 2 a it follows that 



kl{p,V,v') = « g ln 



(111.25) 



Non-relativistic WDM 



From the expression of the conformal energy (1111.120 and the physical velocity V(p, r) 
in (1III.20j) we see that the particle is relativistic if p ^> m a(rj) and non-relativistic for p <C 
ma(rj). Since the comoving momentum is integrated over and weighted by the distribution 



function, we define 



a NR 



( P 2 )h 



m a 



(111.26) 



eq 



where the average is taken with the distribution fo(p) as the value of a that determines the 
transition between the relativistic and non-relativistic regime, the particle is relativistic for 
a <^a NR and non-relativistic for a > a^R- When the particle is non-relativistic 



V(p, 7]) 



P 



ma{rj) 



therefore 

Writing p = yT od we find 



a N R = (V 2 (t eq ))* . 



7.59 x W~ 4 \/y 2 



C=^keVx /2\i 
m J \g d y 



(111.27) 
(111.28) 

(111.29) 



A weakly interacting massive particle (WIMP) (CDM ) of mass ~ 100 GeV and T d ~ 
lOMeV features (V 2 (t eq )}2 ~ 4 x 10~ 8 , whereas for a WDM candidate with m ~ keV 



1 The condition z <C sj y 2 a is equivalent to (p 2 ) ^> m 2 a 2 (i]), where the average is with fo(p). 
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we find {V 2 {t eq ))z < 10 3 , namely all these DM candidates are non-relativistic at t eq with 
(V 2 (t eq )} < 1. Since the WDM particle is non-relativistic at the epoch of matter-radiation 
equality cl^r 1, we find from eqns. (111.81111.101) that 



Vnr 



V2 



k 



(V 2 (t eq )) 



(111.30) 



"i 



for r] 3> riNR the particle is non-relativistic and relativistic for r\ <C rj^R- 

Since in the non-relativistic stage the physical velocity is given by flIII.27p . the integral 
in (lIII.20j) is easily performed by changing integration variable from rj — » a, we find 

kl(p,T),T}') = y a[u — u'] , (III. 31) 

where we introduced 



u(rj) 



1 



In 



-5- 


V 




4?7* + 77 



u N r < u{rj) < 



(111.32) 



where a NR = a(r] NR ), normalized uirf) so that u(oo) = and introduced 

u NR = In 



V a NR 



(111.33) 



During the radiation era when the WDM particle is non-relativistic, a<C 1 we find that 

(111.34) 



kl(p,ij,ij') = In 



1 

IT]' 



which reproduces the result ( Till. 251) . During the matter dominated era for a 3> 1 it follows 
that 

1 277* 

u(ri) j== — . (111.35) 



^(77) 



77 



Free-streaming wavevector from fluid analogy 

In analogy with the Jean's wavevector in the fluid description of perturbations during 
matter domination, we introduce the comoving free-streaming wavevector 



where 



* (V 2 (t)) 



Pm(t) = ; {nt)) - <*<°» 



a?(t) 



a 2 (t) 



and the value of the velocity dispersion today is 

m 



We note that 



k 



eq 



2tt ^3 

" T" (y 2 {t eq ))i 
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(111.36) 



(111.37) 



(111.38) 



(111.39) 



Therefore for these particles 

k fs {a eq ) > k eq . (111.40) 
We define the free streaming wavevector as 

11.17 ( m \ fg d \ I ! 



A* = k fs (a eq ) = — {—)\^y (Mpc) -1 . (HI.41) 

yy 2 

This scale will be seen to play a fundamental role in the (DM) transfer function. 

For a m ~ keV sterile neutrino produced non-resonant ly by boson decay (BD) that 



decoupled near the electroweak scale[57j] ~ 100), it follows that 

£f s D ~ 14.12 (Mpc)" 1 , (111.42) 

whereas for a similar mass sterile neutrino produced non-resonantly via the (DW) mechanism 
near the QCD scale (ga ~ 30) we find 

kf s w ~ 7.7 (Mpc)" 1 . (111.43) 

and for a WIMP of m ~ lOGeV that decoupled thermally at Td ~ lOMeV one finds 
kf s ~ 10 6 (Mpc) -1 . We will see later that kf s determines the scale of suppression of the 
transfer function. 

It is convenient to introduce 

l^a= ^ = v / 6% = 2V2^(f 2 (W)^. (111.44) 

Kfs A K e q 

where y 2 is given by flIII.23j) for the DM species considered here, and Aj s = 2n/kf s . The 
dimensionless ratio k will be important in the discussion of non-relativistic DM. 
From (llll.3imil.32p we find 



fc/(p,77o,^) ^ya\n[V2 + l] (111.45) 

where 7] is the conformal time today, namely l{p,r)o,rj eq ) is the free streaming distance 
traveled by the non-relativistic WDM particle from matter-radiation equality until today. 
Combining this result with eqn. (1III.39I) we find 2 

l(p, no, V e q ) ^ 0.344 -^=X fs . (HI.46) 



From which it follows that during matter domination A/ s , which is the equivalent of the Jeans 
length for collisionless matter perturbations, is simply related to the free streaming distance 
traveled by the non-relativistic particle moving with average comoving momentum a/ (p 2 ) 
from the time of matter-radiation equality until today, namely A/ s ~ 2.9 (p 2 ),T] , rj eq ). 



2 The slight discrepancy with the result in ref. [56j can be traced back to the difference between matter only 
and matter-radiation evolution. 
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From (1111.290 and (IIII.33D we find 



u NR ~ -4.27+ - In 



keV 



50\i 



(111.47) 



where the argument of the logarithm is 0(1) for m 
the (DW) or (BD) mechanisms. 

From eqns. fill. 8|IIII.30j) we find the relation 



keV sterile neutrinos produced via 



Vnr 



V2 



> eq 



{v\t eq )y. 



V2jb 



hence from the definition of k, eqn. (1III.44|) and (1III.30|) it follows that 

k Vnr = - ■ 



(111.48) 



(111.49) 



Therefore comoving modes that entered the horizon when the particle is still relativistic 
correspond ton>2^>k>kf s whereas those that entered when the particle is non- 
relativistic correspond to k < 2 ^> < fcj s . The main corollary is that the free streaming 
wavelength is of the order of the size of the horizon at the time when the (DM) particle 
transitions from being relativistic to non-relativistic. 

This is important: when the particle is relativistic the free streaming distance grows with 
the comoving horizon 77 and free streaming is most efficient to erase density perturbations, 
whereas when the particle is non-relativistic, the free streaming distance grows only with the 
logarithm of the comoving horizon and free streaming is less efficient to erase perturbations 
because the particle free streams with a small velocity. Therefore the dimensionless ratio k 
indicates the regimes in which free streaming is more (k 2) or less (k <C 2) efficient to 
suppress density perturbations. 



B. Initial conditions 

Initial conditions are determined deep in the radiation dominated era and when the 
wavelengths are well outside the horizon. We will only consider adiabatic initial conditions 
for which all the radiation components feature the same Sp r /p r and (non-relativistic) matter 
perturbations obey 




(111.50) 



For the radiation component temperature perturbations correspond to a perturbation in the 
distribution function 

F 1>r (^r h ) = -Q(k i r H )p^^-) ; 8(fe,^) = (111.51) 

so that 

=4e(£,7ji). (111.52) 

/ i,r 
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For superhorizon perturbations when perturbations in the radiation component are nearly 
constant the temperature anisotropy is determined by the Newtonian potential (71-731 

e(k,T )i ) = -- ( j> i (k) 5*77* «1. (111.53) 
Initial conditions for adiabatic perturbations of the matter component also correspond to 

'dfo,m(p)\ 



F hm {kJ ]m ) = -e(£»fc)p(-^W) (IH.54) 



which leads to 



«M<) _ _ e( M,, J f(f^A* -3e(M,) = I (>-> , ,iii,5, 



Pm / P 2 fo,m(p)dp ' 4 \ p 



t.r 



The subtlety for WDM candidates is that in setting up initial conditions for superhorizon 
fluctuations, small comoving scales are superhorizon when the WDM candidate is rela- 
tivistic and intermediate and large comoving scales are superhorizon when the particle has 
become non-relativistic. However, adiabatic initial conditions for all modes are determined 
by flIII.54j) . Indeed, when the WDM candidate is relativistic such initial condition yields an 
energy density perturbation which is adiabatic for a radiation component and when the par- 
ticle is non-relativistic it gives the corresponding relation ( 1III.50I) . Therefore adiabatic initial 
conditions for all modes (superhorizon at the initial time rji) for the WDM perturbations 
are 

F 1 (k,p; Vl ) = ^ l (k)p(^-) ; *ifc<l, (IH.56) 

where fo(p) is the unperturbed distribution function for the DM candidate, and <f>i(k) is the 
primordial gravitational potential determined during inflation. 
In what follows it is convenient to define 

n(t -> \ Fi&PiV) 7( \ fo(p) 

F{k,p;ij) = ; f(p) = (111.57) 

n n 

where 

no = J ^5/o(p), ( m -58) 
is the density of (DM) today. Furthermore, we introduce 

= f -0j3 p (k,r]). (ni.59) 

which becomes 5p m / p m after the DM particle becomes non-relativistic, its initial condition 
is 

6i(k) = 6(k,T)i) = ~^>i{k) ; for k Vi < 1 . (111.60) 
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C. Long wavelength perturbations: 



We begin by studying the evolution of 4>(k, rf) for long-wavelength modes that remain 
superhorizon throughout, to establish the normalization of the transfer function. 

For k — >■ the solution of the Boltzmann equation (1111.191) becomes the same for DM or 
radiation (relativistic) components namely 



F( V ) = F( m ) -[p^\ [<Kv) - <!>{*)] 



(111.61) 



where we have suppressed the argument k since we consider only k = here. 
For the radiation component we write, following eqn. flIV.2p 



(111.62) 



leading to the solution 



Q(ri) = </>(rj) - -<j>. 



(111.63) 



where we used the initial condition ( 1III.53I) . For DM perturbations, from eqn. (1III.59I) we 
obtain 

5(r)) = 3<P( V ) - |& (111.64) 

where we used the initial condition ( 1III.60I) . 

For a DM particle that decouples while relativistic and during the stage when it is still 
relativistic 5p/p ^ 5. However, for a WDM particle with m ~ keV it follows that Sp/p = 5 
for a > o>nr ~ 10 -3 . Hence, for a > a^R the Einstein equation flIII.17p becomes 



2k 2 ~a 2 



3 k 2 q 



+ (1 + 51 



--[55 + 49] 



(111.65) 



where we have used (IIII.14I) . Using the solutions of the Boltzmann equations (IIII.63|IIII.64j) 
for k = 0, and defining = 



we find 



+ 



5 a + 6 


3 


"35 + 4" 


.2 5(1 + 5)_ 


4a 


_ 1 + a _ 



the solution of this equation is 

0(5) 



y/1 + 5 



o 3 



a 3 
% 



3g/ + 4 1 y 3 dy 

L i + y \ V^+y 



+ c 



a 3 



(111.66) 



(111.67) 



and C is determined by giving 0(5jva)- Since a^R < 10 3 for the DM candidates studied 
here, we will take a^R — > whence (piflNR — > 0) = 1, namely we are assuming that the 
DM particle becomes non-relativistic when the Newtonian potential still has the primordial 
superhorizon value. With this initial condition we find 



10 a 3 



16V1 + 5 + 9a 3 + 2a 2 - 81 
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(111.68) 
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a result that agrees with those found in refs. 72|, |77j. For a -C 1 it follows that 0(a) = 
1 — a/10 + 0(d 2 ) therefore the approximation 0(5^) ~ (f)(0) = 0; is very reliable. 0(a) 
decreases monotonically from 0(0) = 1 to 0(oo) = 9/10, and at matter-radiation equality 
0(1) = 0.963. 

For 7^ the transfer function for the Newtonian potential is defined as 



fea»l) = ^T(fc) ; T(0) = 1 . 



(111.69) 



Whereas long wavelength perturbations in the gravitational potential remain nearly con- 
stant, short wavelength perturbations fall off as a consequence of suppression by free stream- 
ing. 

For ka 3> k eq the first term in the left hand side of Einstein's equation ( 1III.17P dominates, 
leading to Poisson's equation 



3 k 2 
4>(k, a) — — 



4 k 2 a 2 



5p 



+ 



6_p 
P 



(III. 70) 



IV. EVOLUTION OF DENSITY PERTURBATION DURING RADIATION DOM- 
INATION. 



Although the Newtonian potential is determined by Einstein's equation ( 1111.131) where 
the right hand side also has a contribution from the DM perturbations during the stage 
when they are relativistic, such contribution is negligible because of the perturbatively small 
effective number of degrees of freedom (/3, A ~ 10~ 2 ) as discussed above. 

Hence, during the (RD) era a « 1 the DM perturbations can be neglected, and the 
evolution of the perturbations is completely determined by the evolution of the radiation 
fluid. In this case there is an exact solution for the Newtonian potential (72 -77 



(z) = -3<j>i(k) 



,73) cos ^) " sin ^ 



V3 



z \3 
V3 ) 



; z = kr] 



(IV.l) 



where 0j the primordial value of the Newtonian potential determined during inflation. The 
solution fllV.ip reflects the acoustic oscillations of the radiation fluid with speed of sound 
c s = 1/V3. 



A. Relativistic DM: stage I 

During the(RD) stage in which the DM particle is still relativistic, namely for kr] ^/y 2 a 

the free streaming distance l(p, rj, 7/) = 77 — 7/ and v (p, rf) = 1, the integrand in 19[) does 
not depend on p. In this case it proves convenient to write 

F(k,p V ) = -6(M; V)p{^) , (IV.2) 

and we find 

z = kr]. (IV.3) 



Q(k,fj,; 77) = —<f)(z) + e~ 



Ifl z 



dz , ( d<t>{z') 
V dz 1 



lflZ 
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Expanding Q(k,/i; 77) in Legendre polynomials, 



Q(k, M! rj) = ( 2/ + !) v)Vi(fi) (W A) 



1=0 

we obtain 



e,(*; 77) = + i 0i(A;) + 2 jf cfc> (^^).wC* - *0 . ( IV - 5 ) 

where we have taken fc^ = 0. The last term describes an ISW contribution akin to that 
in the temperature perturbations of photons (72l|. We note that if the mode remains outside 
the horizon all throughout the evolution during the (RD) stage in which the DM particle is 
relativistic, namely kr] — z <C 1, it follows that 

Qi(z) = ~<j>i(k)8i,o + O(z). (IV.6) 



The WDM density perturbation 



6{k; v) = \f d V f F &P; V)^f, (IV.7) 



therefore during the RD era when the DM perturbation is relativistic 

S(k; V ) = 3e (z) ; S(k;r H ) = ~<l> i {k). (IV.8) 

The monopole Qq(z) begins to grow when it enters the horizon as a consequence of the 
ISW contribution, it reaches a maximun and damps out as a consequence of (relativistic) free 
streaming. This is understood from the following argument: at early time the derivative 
of the Newtonian potential is negative and its modulus increases, reaching a maximum 
approximately at the sound horizon krj ~ v3 ir, whereas the free streaming function jo(z-s) 
is approximately constant for z ~ s, therefore the integrand receives the largest contribution 
near the upper limit, and the total integral peaks near the sound horizon. However, at later 
times the integrand is strongly suppressed by free-streaming since d<p/ds peaks near the 
sound horizon, but for z tc\/3 the free-streaming function suppresses the integrand. Fig. 

© displays e (z)/e (o). 

Although an analytic expression for the integrals in (1IV.5j) is not readily available, we 



can obtain a reliable asymptotic expansion for z > 1. For this purpose it is convenient to 
integrate by parts the derivative of the Newtonian potential, for z > 1 the contributions 
near the upper limit of the integral s ~ z vanish rapidly and the integral is dominated by 
the small s region since the Newtonian potential oc 1/s 2 for large s. Using the asymptotic 
expansion 



M*)= z 2 (iv.9) 

and setting z — > 00 in the upper limit of the integrals we find for z ^> 1 



e t (z) ^ 3<j> t (k) 



siu [z — y] 




+ - , (IV.10) 
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e„(z) 
e„(0) ™ 




z 9 (z) 

e„(o> 




FIG. 1: Left panel q°|q| , right panel: zQq(z)/Qq(0) compared to the asymptotic form pv.iop for 
the monopole. 



this damped oscillatory behavior emerges for z > 15. 

We note that the oscillations in (HV.lOp do not feature the frequency corresponding to 
sound waves, the only remnant of the acoustic oscillations of the radiation fluid in the 
asymptotic form is in the terms featuring the a/3 in the prefactor of the asymptotic form 

poo]) . 

An important conclusion of this section is that during the stage in which the DM particle 
is relativistic density perturbations do not depend on the unperturbed distribution function 
and particle statistics. 

When the particle becomes non-relativistic, for modes k 7]nr ^> 1 the asymptotic form is 
still valid, and the monopole features oscillatory behavior 

sin — 

e (fc) oc — 2 -. (iv.ii) 

K 

This oscillatory behavior is a consequence of the acoustic oscillations of the radiation fluid, 
numerically we find that oscillations arise for k/2 > 15 (see fig. ([I])). 



B. Non-relativistic DM: stages II and III 

When the DM particle becomes non-relativistic (NR) e(p,r]) ~ ma(i]) ; v(p,r)) 
p/ma(r]). It proves convenient to change from 77 to a new variable s defined by 



ds = —P~r =5- s(ri) 



2u(rj) 2V2 



11 



2 k e qd e q 



(IV.12) 



where u(t]) is given by eqn. (1III.32j) . The solution of the Boltzmann equation for the 
normalized perturbation ( 1III.57j) is 



F(k,p;s) = —<p(k,s)(^p ( ^- S j+ \ ds \ ima 2 {s) (j){k\ s')[k -V p f 



snr. 



1 + 



P 



m 2 a 2 (s') 



+ e 



df 

F(k, p; r] NR ) + 4>(k, r) NR ) ( p — 



(IV.13) 
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The initial "time" snr = s(t]nr) corresponds to the (conformal) time at which the DM 
particle becomes non-relativistic. For WIMP's that decoupled thermally for Td <C m at 
conformal time r/d ~ 10 pc during the (RD) era, snr can be taken to be snr = s(Vd)- 
Modes with comoving scales much larger than rjd where outside the horizon at s^r, for 
these modes the initial condition is given by eqn. 56j) . namely 

F(k,p;r] NR ) = \Hk)p (^) • (IV.14) 

On the other hand, WDM particles with m ~ keV WDM decouple when they are still 
relativistic, namely ^> m. For these candidates comoving scales that enter the horizon 
during the (RD) stage when the WDM particle is still relativistic evolve until the particle 
becomes non-relativistic at r\ = t]nr as described in the previous section. Therefore snr = 
s(vnr) and 

F(k,P] Vnr) = -Q(k,K Vnr)p{^-) , (IV.15) 
where Q(k, //; tinr) is given by equations f ]IV.4|IV.5j) with 77 = rj^R. Integrating eqn. (1IV.13j) 



by parts in s' and p, and neglecting the term (p/ma(s)) 2 <C 1 in the non-relativistic limit, 
the evolution of the density perturbation is given by 

6(k, s) = 30(Jfe, s) - k 2 [ ds'a 2 (s')<P(k } s') (s - s') K(k, s - s') 



K(k, 3-s')= I e'^-^fip) (IV.17) 



+ 

where 

K ^ S ~ S ') = j (27r) 3 

determines the suppression by non-relativistic free streaming and 

Q [-ZNR 

S[k,p;r} NR ] = -(j> i {k)e- i, * z » R + 2in / dz'<f){z') e^ 2 ™" 2 ') ; z = krj (IV.18) 

is the result of evolution during stage I and determines the initial condition for the evolution 
during the non-relativistic stages II and III. 

Since f only depends on p, using eqns. flIV.12|IIII.22j) it follows that 



K(k,s-s')=K[a(u-u')] = ^ J y 2 f (y) Jo [ya(u-u')]dy ; N = J y 2 f (y)dy (IV.19) 

and j is the spherical Bessel function. 

The first line in ( 1IV.16j) integrates the gravitational potential during the stages in which 



the particle is non-relativistic. As described above, there are two distinct epochs: when the 
gravitational potential is dominated by perturbations in the radiation fluid and when it is 
dominated by dark matter perturbations. The crossover between the two stages occurs at 
a scale s* = s(a*) that is determined self-consistently, for s > s* the matter perturbation 
dominates the gravitational potential. 

It is convenient to separate the contributions to the gravitational potential from the DM 
and radiation components, writing in obvious notation <p(k,r]) = (f> r (k,rj) + 4> m (k,r)) where 
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4> r (kr]) is given by (llV.lj) . The contribution from DM is obtained from Einstein's equation 
( 1III.17I) which for a > a* reduces to the Poisson's equation for all scales smaller than a few 
Mpc, namely 

0^,10 = -5 Jk (IV.20) 

4 k z a 

For s > s* the integral in flIV.16j) can be split up into the integral from s^r up to s* which 
is dominated by (fi r and corresponds to stage II, and the integral from s* up to s in which 
the gravitational potential is dominated by the DM component (II V. 201) . 



Therefore for s > s*, the density perturbation 5 obeys Gilbert's equation[56|, |57J, |79|, [80 
S{k, s ) = ~^ 5{k, s) + ^ Hi n m J ds'(s - s')K(k, s - s') a(s') 5{k, s') + I[k, s) , (IV.21) 
where the inhomogeneity 

I[k, s] = 3(p r (k, s) — k 2 / ds'a 2 (s')4> r (k, s') (s — s') K(k, s — s') 

J S NR 

and 4> r is the radiation contribution to the gravitational potential given by (llV.lj) . Thus the 
inhomogeneity incorporates the past history during stages I and II. 



C. Kernels for CDM and WDM: 

The kernel K (k, s — s') determines the suppression of WDM perturbations by non- 
relativistic free streaming and depends on the distribution function f(p). For WIMPs 
(CDM) / is the Maxwell-Boltzmann distribution function given by eqn. (IIII.9j) whereas 
for (DW) or (BD) WDM particles f(y) is given by eqn. (IIII.7I) or (1III.8I) respectively. 

1. CDM: Maxwell-Boltzmann distribution function 
For CDM we find 

K(k, s-s') = e~^ u - u ' )2 , (IV.23) 

where u(rj) is defined by eqn. ( IIII.32p . and from the definitions (lIII.44|III.22p . along with 
eqn. (1111.230 . we find 



V6k 

K 



k 



fs 



n on , /100 GeVy /10 MeV\^ / 2 V . . , nrn ^ 



2. WDM: DW distribution function 



With the distribution function ( 1111.70 one finds [56|, |57J, [8 



4 00 (—-\\n+l n 

K { k,.-S) = m = m j:L±— (IV.25) 
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where 



Q = a (u — u') 



0.68 k 
a = — = 0.278 k 

k fs 



(IV.26) 



3. WDM: BD distribution function 



With the distribution function ( 1III.8j) one finds 57 

V2 



K(k,s-s') = K[Q] 



where in this case 



v^Cls) n=1 ( p n ^ 2 





2 


2n + p 


l + n - 




p 




n + p 



Q = a{u — v!) ; a 



0.84 A; 



k 



0.343 k 



fs 



Vn 2 + Q 2 , (IV.27) 



(IV.28) 



We note that in all the cases considered here, the kernels K are functions of the combi- 
nation K 2 (u — U 1 ) 2 . 

The free streaming kernels are suppressed, either exponentially (MB) or as high inverse 
powers (DW,BD) of the ratio k 2 /k 2 s . 



V. COLD DARK MATTER 

For a WIMP of m ~ 100 GeV decoupling at ~ 10 MeV (for which ~ 10) comoving 
scales A 3> % ~ 10 pc entered the horizon well after decoupling and when the particle is 
non-relativistic, in which case we can set rj nr ~ and 



Q{k,p; rj NR ) = -<pi(k) . 



(V.l) 



For these CDM particles, A/ s < 1 pc and for comoving wavelengths A ^ 10 pc it follows that 
ft <C 1 therefore K ~ 1, this amounts to setting (V 2 q )^ = 0, consistently with CDM. The 
perturbation equation ( 1IV.16j) simplifies to 



5(k, s) = 3(j)(k, s)-k 2 f ds'a 2 (s')<P(k, s') (s -s')-- ^(k) 

This equation can be recognized by taking d 2 /ds 2 of both sides, 

d 2 



(V.2) 



ds 2 



5{k, s) — 3(f)(k, s) = —k a (j)(k, s) 



(V.3) 



using d/ds = ad/drj and a/a = 1/r] during (RD) we find 

5 



6 + 



7/ 



+ -4>- k 2 (, 
V 



(V.4) 



which is the equation obeyed by CDM perturbations during the (RD) era [72 
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During this era when is determined by the radiation fluid a 2 (77) = [HlVL m a eq ]rf ; s{rf) - 
\n(rj) /[HlVL m a eq ]z + constant and (p(k,rj) is given by eqn. (1IV.1I) . and eqn. (1V.2I) becomes 



<y(fc,7?) = 9&(*0' 



x cos x — sin x 



X' 



+ 



da: In — -r— 



X NR 



x' ) dx' 



x 



2 



(V.5) 



where x = kr]/y/3. For WIMPs and perturbations with comoving scales A ^> lOpc we can 
set xnr = 0, leading to the result 



x cost x) — sin (x) 



+ - I - Ci(x) + ln(x) + lE } ( V.(>) 

x 2 



where 7^ = 0.577216 • • • and Ci(x) is the cosine-integral function. Fig. ([2} displays 
6(x)/8(0) vs. x = kr)/\/3, where 5(0) = —3(f>i(k)/2. The density perturbation receives 
a "kick" upon entering the horizon at kr] ~ 1. We find numerically that 



<5(x) 
5(0) 



6 ln(x) + 7s ] for x > 10 . 



(V.7) 




We can now estimate the crossover scale at which the Newtonian potential is determined 
by radiation or CDM perturbations. For a C 1 deep in the RD dominated era and for 
subhorizon modes krj ^> 1 Einstein's equation ( 1111131) determines that 

y j ~ 6&(Jfe) cos(x) (V.8) 

Taking the asymptotic behavior (IV. 71) for 5, the Newtonian potential determined by Ein- 
stein's equation (1III.13I) begins to be dominated by matter density perturbations when 

|aln(x)>l. (V.9) 
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For comoving scales smaller than a few Mpc we find that the crossover scale from radiation 
to matter perturbations dominating the gravitational potential is 



a* < 0.1 



(V.10) 



During RD, x ~ d\/2 k/y/3 k eq , therefore for all comoving scales smaller than a few Mpc 
the crossover to the domination of the Newtonian potential by DM density perturbations 
occurs within the RD dominated era. 

Passing to the variable u defined by eqn. (1IV.12j) . for u > u* Gilbert's eqn. (1IV.2ip now 
becomes 

5(k, u ) = ~^5(k,u) + 6 / du'{u-u')a(u')5{k,u')+I[k,u] (V.ll) 



4 k 2 a 



where 



I[k, u] = 3(f) r (k, u) 



8 k 2 [ u * „ 9 
— - \ du'a 2 (u')(f) r (k,u') (u - u') - -(/)i(k) . 
k eq Ju NR * 



(V.12) 



For k 3> k eq and ka 3> k eq which is valid for modes well inside the horizon when DM 



density perturbations dominate, we can safely neglect the first term (IV. lip and because 
during radiation domination kr\ = y/2ka/k eg and for modes deep inside the horizon <p r ~ 
cos(kr]) j k 2 rj 2 we can also neglect the 3<f) r in I[k,u]. We then notice that I[k,u] is linear in 
u and ( IV. lip can be turned into an ordinary homogenous differential equation, 



d 2 
du 



S(k, u) — 6 a(u)8(k, u) = , 



with the initial conditions 



5(k,u*) = I[k,u*] 



d8(k, u) 



du 



dl[k, u] 



du 



Since the variable u depends solely on the combination 

1 



C = v 71 + W) 



tanh[— u) 



(V.13) 



(V.14) 



(V.15) 



(see eqn. 321) ) it proves convenient to write the differential equation ( IV. 13ft in terms of 
C. We find 

A L _ c 2 ) — 
d([ y ^ 1 d(_ 



+ 68 = 0. 



(V.16) 



This is Legendre's equation of index v = 2 with solutions 



A(C) = ~ (3C 2 -i) 

Q 2 (C) = ^(3C 2 -l)ln 



C + i 

c-i 



(V.17) 
(V.18) 



In terms of a rather than ( eqn. (IV. 16ft becomes 

d 2 5 (2 + 35) d5 3 5 

7~ O 



da 2 25(1 + a) da 2 a(l + a) 



(V.19) 
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this is Meszaros' equation [65N67]. We find remarkable that in terms of the variable ( 
Meszaros' equation is simply Legendre's equation of index v = 2. 
The general solution is 



S(k,a) = S g (k)P 2 (C) + 5 d (k)Q 2 (0 ; C = v / T 



(V.20) 



The coefficients 8 g ^ must be obtained from the initial conditions (jV.14|) and the Wronskian 
of the independent solutions P 2 , Q 2 - However, we recognize that the asymptotic solution 
(IV. 7p can be written as 



5(k, a) ~ 65i 



In 



V2k e 7E -i 



V3k 



+ In 



c 2 -i 



(V.21) 



where we used the relation 77 = \/2a/k eq valid during the RD dominated era for 77 <C r) eq 
corresponding to o < 1. Matching (1V.20|) to ( jV.2ip for ( ~ 1 we find 



5 d (k) = -12Si(k) ; 5 g {k) = 68 i {k) In 



V3k 



eq 



(V.22) 



For a> 1 the growing solution is given by S g P 2 ((), namely 



5{k,a) ~ 95i{k) In 



4\/2fc 



(V.23) 



and the gravitational potential becomes for a ^> 1 

9 

<P{k) = —<Pi{k)T CDM {k) , 
where including the long-wavelength normalization (lIII.69j) we find 

TcDM\ 



(V.24) 



OA 45 ^1 



4 a/2 A; e 7E ~ 2 



a/3 fe 



eq 



(V.25) 



is the CDM transfer function for k 3> k eq . This result agrees with that of Weinberg|8l| and 
Wu and Sugiyama 82J and numerically agrees to within few percent with the numerical fit 
provided by Bardeen eiaZ.[83| for k ^> k eq . 

An alternative derivation of this result which is relevant for comparison with WDM below 
begins by defining a new variable 



obeying 



with initial conditions 



d* 
du 



A(k,u) = 5(k,u) - I[k,u] 



A(k, u) — 6 a(u)A(k, u) = 6 a(u)I[k, u) 



A(k,u* 



dA(k,u) 
du 



(V.26) 
(V.27) 
(V.28) 
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Therefore from the solution of (1V.27IIV.28P we find 



8{k, u) — I[k, u] + 6 / a(u')I[k,u']g(u,u')du' 

J u* 



where 



G(u,u') 



W 



P{u)Q{v!) - P{u')Q{u) 



and G[u, u'] = Q(u, u')Q{u — u') is the retarded Green's function obeying 

d 2 



du 2 



6 a(u) 



G[u, u'] = 5{u - v!) 



(V.29) 



(V.30) 



(V.31) 



The functions P{u) = P2(C( M )); Q( u ) — Q2{({u)) are the growing and decaying homogeneous 
solutions of 



du 2 



6 d(u) 



P(u) 
Q(u) 







(V.32) 



and W = 1 their Wronskian. It is straightforward to prove that the solution flV.29j) is exactly 
the same as flV.20j) after using the homogeneous differential equation flV.32j) for P 2 , Q2 and 
twice integrating by parts in u'. 

Since the source I[k, u] remains bound as u — > _ (a — > 00), it follows that asymptotically 
for a 3> 1 

6 r° r° 

5(k,u)^—P(u) Q{u')a{u')I[k,u']du ' = 9a(u) / Q 2 (u) a(u) I[k, u]du . (V.33) 



From fllV.20p and ( 1111.690 we find 



30 k 2 
TcDM^k) = — — 



4 k 2 Uk) Ju 



Q 2 (u')a(u')I[k,u']du' 



(V.34) 



The main reason for describing this alternative in detail is because the form (IV.34j) 
generalizes to the WDM case. 



VI. WARM DARK MATTER: 



Passing to the dimensionless variable u in (IIV. 16j) . eqns. (IIV.21|IV.22p become 

8k 2 ru 



5(k, u) = 3<J)(k, u) — 



rv k 

,v eq J upfn 



a 2 (u')(j)(k,u')U[a(u-u')] du' + 



1 r 3, fdfo(y)\ J 3 . r 1 



N 





znr 



dy 



+ 2 / dz'(f)(z')j 1 [y a(u - u NR ) + z NR - z'] 
Jo 



(VI.l) 
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where 

1 f°° 

U[a(u - u')] = — J y f (y) sin [y a (u - u')] dy = a(u-u')K(k,u-u') , (VI. 2) 

and N is denned in eqn. (IIV.19j) . 

When the DM perturbations dominate the gravitational potential for u > u* which is 
determined self-consistently as explained above, S obeys Gilbert's equation in the form 

5(k,u)-- I a(u')6(k,u')U[a(u- «')] du' = I[k;a;u] (VI.3) 

where we neglected terms proportional to k 2 q /k 2 , and 



I[k;a;u] = 3(j) r (k,u) — 



8k 2 

rv k 

" n eq J u NR 



i r°° 
nJ 



3 , (df (y)' 



a 2 {u') 4> r (k, u') II \a{u — u')] du' + 

3 K 
2 <Pi( k ) jo [y a ( u - u nr) + 2"] 



dz'4> r {z')j 1 [ya(u- u NR ) + ^ - z'] \ , 



(VI.4) 



where we have used z^r = krj NR = k/2. For k ^> k eq the term 3</> r in the first line in (IVI.4P 
is subleading as compared to the second term and will also be neglected in our analysis. 
It is clear from the integral equation (IVI.3j) that S obeys the initial conditions 



8(k, u*) = I[k] a; u* 



d8(k, u) 



du 



dl[k; a; u] 



du 



(VI.5) 



In the first line in (IVI.4I) the kernel II determines the free streaming of WDM perturba- 
tions during the (RD) stage during which the particle is non-relativistic, whereas the last two 
lines are the result of free streaming during the stage when the particle is still relativistic. 
In particular the third term in f ]VI.4j) corresponds to the ISW contribution (IIV.5j) (after an 
integration by parts) studied in section (II V Af) which undergoes damping by free streaming 
during the non-relativistic stage. As it will be seen below, this ISW contribution yields an 
enhancement of the transfer function for k < kf s . 

Thus the inhomogeneity I[k;K;u] is completely determined by the past history during 
stages I and II when perturbations in the radiation component dominate the gravitational 
potential. We have made explicit that the inhomogeneity depends both on k and a (or k). 
For fixed wavevector k the CDM limit is obtained by letting m(gd)z — > oo which lets a — > 
(and k — y 0) with fixed k (see the definition (1III.22j) ) and also unr — > — oo (r^i? —> 0). 

At this stage one can proceed to a numerical integration of (1VI.3j) . however in this article 
we will pursue an approximate semi-analytic treatment valid for an arbitrary distribution 
function postponing a full numerical study to a forthcoming article. 

Before studying (lVI.3|VI.4p . we analyze the asymptotic long time behavior as u — > of 
the WDM density perturbation, which is obtained by neglecting the source term / since it 
is bounded in time. 
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For u — y it follows from (1III.35j) that a(u) ~ 1/u 2 . The integrand in (1VI.3j) is dominated 
by the region v! ~ u ~ 0, assuming that 8(k,u) — y 5(k,0)(—u)~ 13 as u — y and using that 
for vl ~ m ~ it follows that U[a(u — u')} ~ a(w — u ') , and we find 

- / a(i*') g(A;, u')U[a(u - u')] du' ~ 8(k, 0) " V ; . (VI.6) 
cy. J u * P\P ~t~ 1) 

therefore there is a self-consistent solution of eqn. ( 1VI.3I) (for 1 = 0) with = 2, —3 corre- 
sponding to the growing and decaying solutions S g (k,u) oc a; 5d(k,u) oc 1/a?/ 2 respectively. 
This is an exact result which shows that asymptotically for a ^> 1 5 oc a. 

The Volterra equation of the second kind ( 1VI.3j) has a solution in terms of the Fredholm- 
Neumann series. However this iterative solution does not make explicit the growth factor 
a exhibited by the exact solution. The analysis of the CDM case in the previous section 
suggests a re-organization of this series that manifestly exhibits the growth factor. For 
this purpose we cast Gilbert's equation ( 1VI.3[) as an integro-differential equation by taking 
derivatives with respect to u. 

The following integro-differential equation is obtained, 

d 2 6 f u _ , , d 2 d 2 

——^S(k,u) — 6a(u)5(k,u) / a(u) 5(k, u) — -U\a(u — u')\ du = — >I[k, u] . (VI. 7) 

du 2 a ./„» du 2 du 2 



Performing the same asymptotic analysis in the limit u — y 0; d(u) ~ 1/u 2 leading to ( 1VI.6I) 
we find in this limit 3 

fi r u (ft i r°° 

--J^a(u')5(k,u')—U[a(u-u')]du'^a 2 y 2 5(k,u) ; y 2 = - J y 4 f (y)dy. (VI.8) 

This leading asymptotic behavior can be incorporated in ( 1VI.7I) by writing 

d 2 — 

[a(u - u')] = -a 2 y 2 U[a(u - u)] + a 2 U [a{u - u)] (VI.9) 

where 

_ If 00 

U[a{u-u')]=- y fo(y)(y 2 -y 2 ) sin [y a (u - u')]dy (VI.10) 
Using the original integral equation f ]VI.3j) we obtain 



d 2 f w ~r 

-j-^5(k, u) — 6d(u)5(k, u) + n 2 5(k, u) — 6 a J d(u) IT [a (u — u')]5( k, u) du 

d 2 

— I[k,u] + K 2 I[k,u] (VI.11) 



were we used the definition (IIII.44j) . 



The last term in the first line in f IVI.llj) can be interpreted as a non-local potential with 



a memory kernel II[q;(m — ?/)] . It is straightforward to show that II[q;(m — ?/)] oc 



u — u') 3 



3 This can be found self-consistently by proposing S(k,u) oc (—u) 13 and following the steps leading to 
(jVLH) 
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as vl — y u and from the results for the kernels (1IV.23|IV.25|IV.27j) that it falls off as a high 
power (or exponential) of the argument for the distribution functions considered here. 

Furthermore, we have already established that asymptotically 5(k,u) oc a oc l/u 2 , imple- 
menting the same analysis leading to (1VI.6P and replacing this asymptotic behavior in the 
memory integral in (1 VI. lip we find that asymptotically as u — y it behaves as 



a[u 



')n[ 



a [u — u 



')] S(k, u) du oc ln(— u) oc ln(c 



(VI.12) 



therefore its contribution is subleading in the asymptotic limit a — > oo as compared to all 
the other terms in the first line of fIVI.lip . 

Hence, we conclude from this analysis that the memory integral in (I VI. lip can be con- 
sidered as a perturbation. 

Again, it is convenient to introduce the combination A(k, u) given by (1V.26|) that satisfies 



d 2 

— A(k, u) - 6a(u)A(k, u) + n 2 A(k, u) = Ga(u)I[k, u] + J[5; u] 



with the initial conditions given by (jV.28|) . where 

J[5\u] = Qa I d(u) Y[[a(u — u')]5(k, u) du . 



(VI. 13) 



(VI. 14) 



The solution of (1VI.13P with the initial conditions (IV. 281) is completely determined by the 
retarded Green's function obeying 



d 2 
du 



— 6a(u) + k 2 



G[u,u'} = 5(u-u'). 



(VI. 15) 



The formal solution of fIVI.lip with initial conditions flVI.5P is 



6{k,u) = I[k,u] + / g{u,u') 

J u* 



6a{u')I[k,u'\ +J[5; u'] 



du' 



where 



Q(u,u') 



1 

W L 



P(u) Q(u') - P{v!) Q(u) 



(VI. 16) 



(VI. 17) 



where P, Q are the linearly independent growing and decaying homogeneous solutions of the 
fluid-like equation 



- 6a(u) + k 2 

du z 



P{k;u) 
Q(k; u) 







(VI. 18) 



and W is their (constant) Wronskian. The formal solution (IVI.16P is again an integral 
equation, however it is a re-summed form of the Fredholm- Neumann solution of (IVI.3P that 
displays the asymptotic growth factor explicitly since asymptotically the growing solution 
of (IVI.18P P(k;u) features the growth factor oc a (see below). 

From the analysis above, we note that the inhomogeneity J is subleading compared to 
the first term dl[k,u] for the following reasons: 
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• At early times u ~ u*, J vanishes as [u — u*) 3 whereas al[k,u] remains finite. 

• Asymptotically at long time (u — > 0; a — > oo) al[k,u] ~ al[k, 0] oc a whereas J7" oc 
ln(a). 

• At long wavelengths k — » for which a — >■ (k — > 0) it follows that J — > 0. This is 
the CDM limit. 

• For short wavelengths free streaming suppresses density perturbations, this is manifest 
in the expression (1VI.4|) . In an iterative solution 5 is suppressed by free streaming and 
the term J involves a further suppression by the kernel II with respect to I. 

Hence the term J can be treated perturbatively as argued above, giving rise to a systematic 
Fredholm-Neumann iterative solution of (1VI.16j) formally in powers of the free streaming 
kernels II which for (WDM) are strongly suppressed by large inverse powers of k at small 
wavelength (see the expressions ( 1IV.23tnV.27j) or exponentially suppressed as for (MB) (see 

flrZ23D) 

5(k, u) = 5 i0) (k } u) + 5 {1 \k, «) + ••• (VI.19) 

where 



S (0 \k,u) = I[k,u] +6 / Q{u ,u')a(u')I[k,u'} (VI. 20) 

J u* 

5W(k,u) = / giu^u^J^-^-u^du' -n > 1 (VI.21) 



note that 5^(k,u) is first order in the free streaming kernels, 5^(k, u) second order, etc. 
We refer to the zeroth-order solution (IVI.20j) as the Born approximation because of its 



similarity to quantum scattering theory. In references [56|, |57| it was shown that the Born 
approximation is reliable in a wide range of scales. In what follows we will study the transfer 
function in the Born approximation as a prelude to a full numerical study of (IVI.3j) and its 
comparison to the Born and higher approximations to be reported elsewhere. 

We note that the Born approximation is exact for CDM since in this case a = (conse- 
quently n = 0). 

It remains to obtain the homogeneous solutions P, Q of the fluid-like equation (1VI.18j) . 
which becomes more familiar when written in terms of cosmic time t, 



d 2 nTT d k 2 (V 2 (t)) , 




(VI.22) 



where p m (t); (V 2 (t)) are the density and the velocity squared velocity dispersion of the 
DM particle given by (1III.37|IIII.38I) . This is equivalent to the Jean's fluid equation for 
non-relativistic matter recognizing that k/a(t) = k p h ys {t) is the physical wavevector, and 
replacing the (adiabatic) speed of sound by the DM particle's velocity dispersion. The term 
proportional to k 2 plays the role of a pressure term and its origin is traced back to the 
free-streaming kernel II in Gilbert's equation (lVI.3p . 

We emphasize that whereas the fluid equation (1VI.22|) suggests acoustic-like oscillations 
and is familiar, it is only half the story, it has no information on the suppression of perturba- 
tions by free streaming. The solution of Gilbert's equation flVI.19|VI.20|VI.21j) is completely 



determined by the inhomogeneity and initial conditions, these are determined by the past 
history and describe the suppression of density perturbations by free-streaming. 
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A. Meszaros' equation for WDM 



Just as in the CDM case (see equations flV.13|IV.16j) ). it is convenient to pass to the 
variable (, in terms of which the homogeneous equation flVI.18p becomes 



d 



d 



i-a-j^-^-TT + ^o + i) - 



d( 



d( 



i-C 2 



; v = 2 



(VI.23) 



this is the associated Legendre equation with indices v = 2 ; in. We choose the growing and 
decaying solutions respectively as 



P( K ,C) = Re< 



c-i 

C + i 



2.3.1 iKj . 



i-C 



(VI.24) 



c) = S inh(™) I _ k) F 



2,3;1 



i-C 
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where F[a, b; c; z] is the hypergeometric function. We find 

P(k, u) = cos(k u) F r (k, C( u )) + K sin(K u) H(k, ((u)) 



(VI.25) 



(VI.26) 



Q(k, u) 
where 



1 



3P{k,u) + (k 2 - 2) 



cos(ku) H(k, C(u)) ^ — F r (k, C(u)) 



K 



F(k 3(1 ~ C) I 3(2 -^-O 2 

F R (K, C(«)) - 1 - + (1 + /?2)(4 + fi;2) 



H(k,C) 



3(1-0 



;i + K 2)(4 + ft . 2) 



1 + k 2 + 3 C 



tanh[— u\ 



(VI.27) 



(VI.28) 



(VI.29) 



It is straightforward to confirm that P(0, Q = f^O j Q(0, C) = $2(0 are the Legendre 
functions solutions of Meszaros's equation flV.17|V.18p for CDM perturbations. In fact, in 
terms of the variable a equation ( 1VI.18P (or alternatively eqn. ( IVI.23P ) becomes Meszaro's 
equation for WDM, 



d 2 (2 + 35) d _ 

da? 2a(l + a) dd 2a(l + a) ' 4a 2 (l + a 



P 

Q 



(VI.30) 



whose growing and decaying solutions are given by (1VI.26fVI.27p respectively. 

The asymptotic behavior of the growing and decaying solutions for a> 1 ; u — > are 



P(k,u) -»■ 
Q(k,u) -»■ 



3(2 -k 2 



M 2 (l + K 2 )(4 + K 2 ) 

-m 3 (1 + k 2 )(4 + k 2 ) 
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(VI.31) 
(VI.32) 
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from which we extract the Wronskian 



W 



2-k 2 



Therefore we find 



G(u,u') 



2-K 2 



Q-MODES 
k = 6.3 (fundamental) 




^ 1 

ST 





k=13.6 



K=10 



/) - P{k,u')Q{k,u) 




P-MODES 




K = 10 






N \k=13.6 


1 K = 




- 





(VI.33) 



(VI.34) 
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FIG. 3: Mode functions of fluid equation (|VI.18p . Q(k, u) are the decaying and P(k, u) the growing 
solutions. The "fundamental" decaying solution features a node at matter-radiation equality. 



For a ^> 1 when the gravitational potential is determined by DM perturbations, using 
Poisson's equation flIV.20j) . the definition of the transfer function flIII.69j) and the solution 
for 5 (IVI.16j) along with the asymptotic behavior flVI.31j) of the growing solution P(k, u) 
leads to an exact expression for the transfer function 



WDM 



(k; k) 



5 kl q 



Q{k,v!) 



6a(u')I[k;K;u'] + J[S;u'} 



du' . 



fc 2 (l + K 2 )(4 + K 2 )0i(A;) J u * 

(VI.35) 

The CDM transfer function Tc DM (k) corresponds to setting a = 0;r] NR — >■ which sets 
K = 0; Unr — > — oo and J = 0. In the Born approximation we obtain 



T B (k; k) 



-30k 2 eq 



fc 2 (l + /€ 2 )(4 + /€ 2 )0,(fc) J u « 



Q(k, u') a(u')I[k; k; u'}du 



(VI.36) 



and as explained above the Born approximation is exact for CDM (for k 3> k eq ). 

TcDhiik) is given by (IV.34j) and its leading behavior for k 3> k eq is given by (IV.25|) . It is 
convenient to normalize the WDM transfer function defining 



T(k) 



WDM 



(k; k) 



TcDM^k) 

where WDM refers to k ^ 0. In the Born approximation we find 

a r r° 

T B (k) 



;i + k 2 )(4 + k 2 ) 



fu* Q( K i u ') o,(u')I[k; k; u'} du' 



J^Q 2 (u')~a(u')I CDM [k;u'} du'_ 



(VI.37) 



(VI.38) 
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where Q2 is the Legendre function given by eqn. (1V.18j) . 

, > 1 



a\u) 



sinh \u] 



and 



tCDM 



[k;u] = l[k;0;u) 



(VI.39) 



(VI. 40) 



The matching scale u* describes the transition from when the gravitational potential is 
dominated by the radiation fluid to when the DM perturbations dominate. In the CDM 
case analyzed in sect ion (jVj) we found that this scale is smaller than the scale of matter- 
radiation equality. From the result (IV. 7j) and the analysis leading to ( 1V.21|) we also found 
that the density perturbation in CDM depends logarithmically on the change of scale and 
for k ^> k eq taking the matching scale 



u 



u 



eg 



hi 



V2 



y/2 + 1 



-0.881 



(VI.41) 



yields a correction which is of order k1 q /k 2 <C 1 in the small scale regime studied here. For 
WDM, free streaming makes the dependence on this scale even weaker, and it is evident 
from the expression (IVI.38P that the contribution from a <C 1 is suppressed. Hence, in 



our analysis we take u* 



u 



cq 



—0.881. A comprehensive numerical analysis confirms the 



insensitivity on the choice of scale for k 3> k eq . 

It is convenient to divide the inhomogeneity (1VI.4|) by — 30j(/c) which cancels in the 
ratio (1VI.36j) . Furthermore since the integrals in (1VI.36j) range from i] eq ^ rj ^ 00 and 
4>(k,r]) oc l/(ki]) 2 we can safely neglect the term 30 r in the first line in (1VL4I) as compared 
to the second term for k 3> k eq . Thus in the ratio (1VI.36I) / simplifies to 
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leading to 

I CDM [k- u] = P ~a\u') (p(k; v!) (u - v!) dv! + | , (VI.48) 

which along with ( IV. 181) determines the denominator in (1VI.36I) . 

In the appendix we provide an explicit form for (1VI.43[) . we gather all the relevant re- 
sults, and provide a concise summary of the Born approximation for an easy numerical 
implement at ion . 

The contribution Iisw is a result of an integration by parts in eqn. flIV.3j) and is the only 
contribution that vanishes in the CDM limit. It originates in stage I during (RD) when the 
WDM particle is still relativistic. 

Figures (T4"f5]) displays the ratio T and its logarithm for both cases of non-resonant sterile 
neutrino production (DW,BD). The production via boson decay at the electroweak scale 
leads to a colder species for two reasons: i) the effective number of degrees of freedom at 
decoupling g d is larger, therefore the particle is colder today and at matter-radiation equality, 
and ii) the distribution function (IIII.8P favors small momenta and yields a smaller velocity 
dispersion (see eqn. (1III.23I) ). This is manifest in the transfer functions displayed in fig. (j4j): 
it is clear from this figure that the wavevector scale of suppression for DW-produced sterile 
neutrinos is smaller than for the BD-production mechanism for the same mass. 




k (Mpc)" 1 k (Mpc)" 1 

FIG. 4: Tg(fc) for DW, and BD for m=l,2 keV. Sterile neutrinos produced via the BD-non-resonant 
mechanism are colder for the same mass. 



B. ISW enhancement: 

As discussed above the contribution Iisw is a direct consequence of the evolution of 
density perturbations during stage I during the (RD) era described by eqn. (1IV.3I) . and 
vanishes in the CDM limit. Therefore it is a distinct contribution to the WDM transfer 
function, and only arises from the time evolution of the Newtonian potential driven by the 
acoustic oscillations of the radiation fluid, i.e. an ISW effect. 

This contribution is "out of phase" with the first two terms Ii^: the Bessel functions 
jo of these two terms are decreasing functions of k until their arguments vanish. Instead, 
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the Bessel function j 1 grows during the initial interval when j decreases. As a result hsw 
grows for small k. This is precisely the behavior displayed in fig. ([1]) corresponding to the 
I = (monopole) component of the density perturbation (1IV.5j) (integrating by parts the 
integral term the jo becomes ji). 

Since the maximum value of r\ during stage I is t/nr and kr/^R = k/2 the analysis following 
eqn. (1IV.8|) suggests that hsw features a peak when the wavelength of the perturbation is 
approximately the sound horizon at r] NR , namely kr] N R ~ y/3n or k w 2n\^3. This analysis 
suggests that hsw features a peak at k < kf s because the argument of the Bessel function 
is now shifted towards the positive values (since u — u^r > 0). The presence of a peak can 
also be gleaned from hsw directly, since for small z' <f(z') is nearly constant but j\ grows, 
featuring a maximum when its argument is ~ 2, which obviously suggests a peak at k « kf s . 
Therefore the hotter species, with smaller kf s must feature a peak at a smaller value of k 
when compared to the colder species which features the peak at a larger value k because 
of a larger value of kf s . This expectation is borne out by fig. ([6]) that displays the ISW 
contribution to the Born ratio T# ( 1VI.38I) . The ISW enhancement extends to larger values 
of k for the colder species for the same mass (BD) as a consequence of a larger value of kf s . 

For small k the contributions I 2 and hsw feature opposite signs, therefore the ISW 
enhancement competes with and is partially cancelled by h yielding an overall suppression 
of the transfer function with respect to CDM. Nevertheless, the ISW enhancement prolongs 
the region in k where the transfer function is closer to that of CDM. 

For k > 30 {k ^> kf s ) the ISW contribution features oscillations as discussed in section 
( 1IV Aj) and shown explicitly in fig. (JTJ). 



C. On the origin of WDM acoustic oscillations: 

The Q and P modes (lVI.27pVI.26j) feature acoustic oscillations as displayed in fig. (j3j), 
and only the Q modes enter in the evaluation of the transfer function (1VI.36j) . This mode 
function always vanishes at u = (today), and there is a particular "fundamental" mode that 
features only one other node at matter-radiation equality, for k ~ 6.3. In the integral leading 
to the transfer function (IVI.36[) the mode function Q multiplies the three contributions to 
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ISW contribution 




I displayed in (IVI.42tiVI.45) . The integral over y with the distribution function leads to the 
dephasing of the oscillatory functions in Ij, I 2 , hsw and their suppression by free streaming. 
However, we can identify some of the more obvious oscillatory contributions. From the 
study in section (1IV Aj) and the results displayed in fig. ([I]), the oscillations from the ISW 
component begin at z NR > 15 (k > 30) or k > 5\/6 kf s , and are suppressed by free streaming 
during stages II) and III). This suppression is encoded in the y integral with the distribution 
function which contributes during the stages when the particle is non-relativistic. 

The explicit form of I\ given in the appendix, ( 1A.7I) reveals at least two contributions that 
lead to oscillations, these are the term sm[ayU]/ay in the first line, and the second line 
in f]A.7|) . After integrating in y these contributions are proportional to the free streaming 
kernels ( 1IV.23IIIV.25IIIV.27I) . however, although these contributions do not feature oscilla- 
tions after the integration in y by themselves, they are multiplied by the mode function Q. 
Therefore the last term in the first line in (IA.7j) leads to oscillations for wavevectors larger 
than that of the "fundamental" Q-mode. The second line in flA.7j) yields a contribution of 
the form 

r sin(x NR )- 

CX 1 

L x NR . 

times a function suppressed by free streaming. With xnr — k/2 this contribution vanishes 
for k <C kf s , reaches the value 1 at n = 2tt and oscillates around one for k ^> 2n. Therefore 
this function reaches its asymptotic value ~ 1 for values of k near the "fundamental" mode. 
This analysis leads us to suggest that oscillations in the transfer function begin when the 
"fundamental" mode is excited, namely k > 6.3. 

For values of k > 6.3 the nodes in the mode functions Q between matter-radiation 
equality and today lead to oscillations in the transfer functions. Therefore we conclude that 
oscillations are manifest for 

k>2k fs . (VI.49) 

This expectation is approximately borne out, for (DW) with kf s ~ 7.7 (Mpc)" 1 we see from 
fig- (El) that oscillations begin at k « 11 (Mpc) -1 and for (BD) with kf s m 14 (Mpc) -1 , fig. 
([8]) shows oscillations beginning at k w 31.5 (Mpc)" 1 . The period of the oscillations is more 
difficult to assess because the various terms are out of phase leading to beating of frequencies 
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(a hint of this is observed in ln(T) displayed in fig. ((7j)), however, the approximate estimate 
k ~ 2k f s for the emergence of oscillations is confirmed by the numerical analysis. 

It is important to recognize that both I\,Iisw originate in the acoustic oscillations of 
the radiation fluid, which couple to the WDM perturbations via the Newtonian potential. 
Therefore in this sense, the origin of the WDM acoustic oscillations at small scales is similar 
to the small scale oscillations in the CDM transfer function obtained in ref.[70|]. In that 
reference the oscillations originated from the direct coupling of the CDM particle to the 
radiation fluid prior to decoupling, whereas in this work the coupling is indirect through the 
gravitational potential and the past history of the evolution during stages I and II. 

At the scale where WDM acoustic oscillations emerge the transfer function is strongly 
suppressed by free-streaming and as a result of this suppression in the power spectrum the 
relevance of these WDM acoustic oscillations for structure formation is not clear. However, 
it is conceivable that the effect of the oscillations will be amplified by non-linear gravitational 
collapse, leading to enhanced peaks and troughs in the matter distribution at low redshift. 

The (comoving) scales for these oscillations k ao ~ 11 (Mpc) -1 for (DW) and k ao ~ 
31.5 (Mpc) -1 for (BD) could lead to dumpiness in the mass distribution with mass scales 
M DW ~ 3 x 10 9 M or M BD ~ 1.8 x 10 8 M respectively. 




k (Mpc)" 1 k (Mpc) 



FIG. 7: Acoustic oscillations at small scales: (DW) species. 



The smaller amplitudes of acoustic oscillations for the (BD) species as compared to the 
(DW) case is consistent with the fact that (BD) sterile neutrinos are colder and feature 
smaller velocity dispersions. 



D. Power spectra: interpolation between large and small scales. 

The power spectra normalized to CDM is given by 

p(k) = \T{k) 

and the full power spectra is therefore, 

P(k) = P CDM (k)P(k) . 



(VI.50) 



(VI.51) 
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FIG. 8: Acoustic oscillations at small scales: (BD) species. 



Since the transfer function for WDM particles is indistinguishable from that of CDM for 
small k, and as we have pointed out above the result (IV. 251) coincides within a few percent 
with the result by Bardeen et. al. fiifor k ^> k eq , we use the numerical fit provided by 
Bardeen et.a/.j83[ for the CDM transfer function (without baryons) to extrapolate PcDM(k) 
to large scales: 
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where A is the overall amplitude and is determined by the power spectrum of scalar fluc- 
tuations during inflation [72|, and n s ~ 0.96 is the index of scalar perturbations during 



inflation |69J. Without baryons and with three relativistic (standard model) neutrinos [83 
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(VI.53) 

Combining eqns. ( 1VI.50llVI.51IIVI.52l) and using the Born approximation for T(k) we find 
the following expression for the power spectra that interpolates between large and small 
scales, 



P{k) = Ak r 
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-0.881 and the mode func- 



The inhomogeneities J, I CDM are given by (IVI.42tiVI.48j) . u 
tions Q2, Q are given by eqns. ( 1V.18PVI.27[) respectively. The appendix gives a simplification 
of these terms along with a numerical implementation. This compact expression provides 
an interpolation between large and small scales that describes accurately the CDM limit 
for long- wavelengths k k/ s and captures the free streaming suppression at small scales 
encoded in the Born approximation. Its numerical implementation is fairly straightforward 
for arbitrary distribution functions, mass and decoupling temperature. 
This is one of our main results. 
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E. Comparison to numerical results from Boltzmann codes: 



The (WDM) power spectrum for non-thermal sterile neutrinos produced via the (DW) 



mechanism has been studied in refs. [3ll.l38l - l41| . The most recent studies using the Boltzmann 



codes CMBFAStQ and or CAMMlJhave been reported in refs.0, 51 E 



The results 

of ref. kl| coincide with those of ref.[3l| and are summarized by the fit given by eqns. (6,7) 
in ref.[3l|. In both refs. 3l|, |41| th e distribution function for sterile neutrinos is that given 
by eqn. (1III.7I) obtained in ref. 371 ] . However, the fitting function eqn. (6,7) given in ref. j3lj 
(which reproduces the results of ref. 4l|) fits the results of the Boltzmann code in the range 
k < 5 h M2C _1 0. 

In ref.|40| the kinetic equation for production of sterile neutrinos given in ref. 37 



was 



solved numerically and the solution was input in the numerical Boltzmann codes. In this 
reference the explicit form of the distribution function is not provided but instead a fitting 
formula for the transfer function normalized to CDM is given, eqn. (11,12) in this reference. 



Whereas both fitting functions in refs. [31l . |40( are of the same form, they differ in the powers 
of momenta: at large k the fitting formula (11) in ref. 40J falls off with a power ~ ^-6.93 
whereas the fit given by eqn. (6) in ref. [31] falls of with a power ~ A; -10 . Therefore at small 
scales there is a large difference between these fits, whereas at large and intermediate scales 
there is a substantial agreement (see fig. 4 in ref.(40j). Because in ref. 40] the distribution 
function has been obtained directly from the numerical integration of the kinetic equation 
derived in ref. 37J , it is not clear whether the main differences with the results of ref. [3l| are 
a consequence of the distribution function obtained numerically and input in the Boltzmann 



code being different from the form T|) which is the one used in refs. 3l|, |41 



Because our study relies on a pre-determined form of the distribution function and we 
neglect baryons, we can most directly compare our results with the distribution function 
( 1III.7j) to the results in ref. 3l| , which also uses the form (1III.7P and neglects baryons, however 
it includes fl\ = 0.7 which our study does not. 

We compare our results for the transfer function T{k) (normalized to CDM) given by 
(1VI.38I) with those obtained from the fit given by eqns. (6,7) (for the non-thermal case) 
in ref. [3jJ, with the caveat that this fit may not be the correct description of the power 
spectrum for k > 5 h Mpc -1 as suggested by the discussion in ref. [HI]. We also compare to 
the fit (11,12) in ref.[40j], although this may not be fair comparison because we assume the 
distribution function (1111.70 whereas in ref. |40j the effective distribution function may be dif- 
ferent and the difference cannot be quantified in absence of a functional form. Furthermore, 
we use the "standard" value = 10.75 for the comparison, whereas as discussed in ref. jloT] 
the actual value may differ because this species of sterile neutrinos is produced very near 
the QCD phase transition where the effective number of relativistic degrees of freedom vary 
rapidly. Recognizing all these caveats we present the comparison of the transfer functions 
normalized to CDM in the range of masses and scales displayed in refs. 3l|, |40| in fig. (Q, 
m = 0.5, 1.0, 1.7 keV: the solid line is T{k) from the Born approximation ( 1VI.38j) . the dashed 
line is the fit given by eqns(6,7) for the non-thermal case in ref. [Hj, the dotted line is the 
fit (11,12) in ref.jloj. 



We find a remarkable agreement, to less than 5% with the fit given by eqns. (6,7) (non- 
thermal case) in ref. 31] in a wide range in which their fit is valid (see discussion in ref. [3l| ) 
for m > 1 keV the agreement is substantially better in a far larger range. In fig. fl9]) the 
comparison is in the range displayed in refs. 3l|,|40| to highlight agreements and discrepancies. 
In all cases reported in the literature the range studied or displayed are for wavectors k far 
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k [h/Mpc] 



FIG. 9: Comparison of the transfer function for DW with = 10.75 normalized to CDM with 
the results from Boltzmann codes. The solid line is the semi-analytic result from eqn. (|VI,38j) , the 
(blue) dashed line is the result from the interpolation eqns.(6,7) (non-thermal case) from ref . ( (3l| ) . 
the (red) dotted line is the result from the interpolating fit eqn. (11,12) in ref. ( 4y]). For all cases 
h = 0.72, n DM h 2 = 0.133,c/ d = 10.75. 



smaller than the range in which the acoustic oscillations become manifest. The approximate 
estimate (IVI49P for the threshold suggests that for m = 0.5, 1.0, 1.7 keV oscillations should 
be manifest for k > 5.4, 10.8, 18.5 (Mpc) -1 (corresponding to k > 7.5, 15.0,25.6 h (Mpc) -1 
respectively). Fig. f[T0~j) displays T{k) from ( 1VI.38I) in a linear-linear scale for k > 2k f s for 
m = 1.0, 1.7 keV. These figures are the continuation of the same T{k) displayed as solid 
lines in fig. ([9]) to the smaller scales k >2kf s in each case. 



DW; m=1.0 keV ; u = 10-75 



DW; m=1.7 keV ; gj = 10.75 




k (Mpc) 



k (Mpc) 



FIG. 10: T(k) from the semianalytic approximation (|VI.38|) displaying the acoustic oscillations 
at small scales k > 2kf s ~ 10.8,18.5 (Mpc) -1 for m = 1.0, 1.7keV respectively. Note that the 
horizontal scale is in (Mpc) -1 and that vertical scales differ by a factor 5 between the two figures. 



This comparison, with all the caveats mentioned above, suggests that the semi-analytic 
formulation along with the Born approximation summarized by flVI.381) captures the essential 
physical processes and provide a reliable tool to study the transfer function and power spectra 
for arbitrary distribution functions. 
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F. Impact on N-body simulations and Lyman-a constraints: 



State of the art N-body simulations of galaxy formation 27], l28| and large high resolution 
data sets of Lyman-a: forest spectraj29-31| have been used to constrain the mass of WDM 



particles [30, 131 



The most recent large scale N-body simulations 27J, |28| incorporate WDM by considering 
a power spectrum that is cutoff at small scales, however, initial velocity dispersion is not 
yet included in the simulations. Extracting constraints from the Lyman-a forest involves 
also large scale numerical simulations, and the most recent constraints [30] on the mass of 
the WDM particle rely either on a thermal or (DW) distribution functions. The (DW) dis- 
tribution function is proportional to a thermal distribution function and the proportionality 
constant only determines the abundance but is irrelevant for the free streaming length or 
indeed the transfer function (as can be gleaned from the previous sections). 

Our study points out that the power spectrum features a quasi- degeneracy in that a more 
massive WDM particle with a (DW) distribution function features a similar power spectrum 
as a less massive one but with a (BD) distribution function in a wide range of scales. To 
make this more explicit, fig. (jlip displays the power spectra normalized to CDM flVI.501) . 
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FIG. 11: The power spectra normalized to CDM for DW and BD with m = l,2keV. Note that 
for the same mass the BD (colder species) is less suppressed than the DW (hotter species). 



From this figure it is clear that P{k) for (DW) with m = 2 keV is almost indistinguishable 
from P(k) for (BD) with m = 1 keV for k < 6 — 8 (Mpc) -1 . This is because the (BD) 
sterile neutrinos are colder for two reasons: they decouple earlier and their distribution 
function favors small momenta, therefore the (BD) WDM particle has smaller velocity 
dispersion. Therefore, we emphasize that the mass is not the only relevant indicator for the 
power spectrum of the WDM particle, but also two important aspects must enter in the 
assessment: the decoupling temperature (the higher, the colder the particle) and the details 
of the distribution function at small momenta: enhanced small momentum behavior leads 
to a colder species and a less suppressed power spectrum, for a given mass. 

Hence the quasi- degeneracy: the current constraints on the mass of the WDM particle, ei- 
ther from (quasi) WDM simulations (quasi because these simulations do not include velocity 
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dispersion of the WDM particle, therefore miss the aspects related to the non-thermal dis- 
tribution functions), or from Lyman-a forest analysis, which typically input thermal WDM 
distribution functions or (DW) distribution function which is indistinguishable from ther- 
mal for the purpose of the transfer function, do not directly apply to non-thermal WDM 
particles. 




FIG. 12: The matter power spectra: P{k) = Ak(T(k)) 2 for n s = 1 (A is the normalization 
amplitude) for CDM, DW and BD for m = l,2keV. Note the quasi degeneracy for DW with 
m = 2keV (d) and BD with m = IkeV (c) in a large range of k < 12 (Mpc) -1 . 



To highlight this point, we obtain the full power spectra for the different species considered 
here using the interpolating eqn. (1VI.54|) . Fig. ( fl2|) displays P(k) for n s = 1 ; fl m h 2 = 0.134 
for the different species considered here. Note how the two cases (c) (BD, m = 1 keV) and 
(d) (DW, m = 2keV) are nearly indistinguishable for k < 6 — 8 (Mpc) -1 . 

Therefore, we conclude that non-thermal distribution functions may evade the constraints 
on the mass of the WDM particles both from current numerical simulations and the Lyman-a 
forest data. 



VII. CONCLUSIONS AND DISCUSSIONS 



In this article we provide a semi-analytic study of small scale aspects of the power spec- 
trum of WDM candidates in a radiation-matter cosmology for arbitrary mass and distribu- 
tion function of the decoupled WDM particle. There are three stages in the evolution of 
density perturbations of WDM candidates that decouple while they are relativistic: stages 
I) and II) describe the evolution during the RD era when the particle is relativistic and 
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non-relativistic respectively but the gravitational potential is dominated by the radiation 
fluid, during stage III, the particle is non-relativistic and matter density perturbations dom- 
inate the gravitational potential. We consider adiabatic initial conditions determined when 
all the cosmologically relevant modes are superhorizon. The collisionless Boltzmann equa- 
tion is solved in the three stages by using the solution at the end of a stage as the initial 
condition for the next stage. The transfer function is characterized by two widely sepa- 
rated scales: k eq ~ 0.01 (Mpc) -1 corresponding to the wavevector that enters the horizon at 
matter-radiation equality and 




where (V^ q )^ is the mean square root velocity dispersion of the WDM particle at matter- 
radiation equality. This latter scale also determines the size of the comoving horizon when 
the WDM particle becomes non-relativistic: 

During stages I) and II) the acoustic oscillations in the radiation fluid dominate the gravi- 
tational potential, leading to an ISW effect that amplifies WDM density perturbations on 
scales larger than the sound horizon at rj nr- This amplification translates in a prolonged 
plateau in the transfer function for k < kf s which is more pronounced for colder species 
since these feature a larger kf s . 

When the particle is non-relativistic and WDM perturbations dominate the gravitational 
potential, the evolution is described by the Boltzmann-Poisson equation which yields an 
integral equation for density perturbations and is equivalent to integro-differential equation 
with an inhomogeneity and initial conditions determined by the past history during stages I 
and II. This equation is amenable to a systematic Fredholm expansion valid at small scales, 
whose leading order is the Born approximation which establishes a direct relation with a 
fluid description of WDM perturbations. The resulting fluid equation is the generalization 
of Meszaros' equation for CDM but with an inhomogeneity and initial conditions that incor- 
porate suppression by free streaming during the first two stages. The Born approximation 
lends itself to a simple numerical implementation for arbitrary distribution functions and 
mass of the decoupled WDM particle. Its main ingredients are the growing and decaying 
solution of the generalized Meszaros fluid equation for WDM perturbations, and the initial 
conditions and inhomogeneity that are completely determined by the past history during the 
first two stages. The solutions of the fluid equations feature (WDM) -acoustic oscillations 
which are manifest in the transfer function and power spectra for k > 2k f s . 

An approximate form of the power spectra that interpolates between large and small 
scales for arbitrary distribution functions is given by eqn. ( 1VI.54|) and a simple and concise 
summary of the main elements of the Born approximation and its numerical implementation 
are provided in the appendix. 

We study in detail and compare the transfer functions and power spectra of sterile neutri- 
nos with mass in the ~ keV range for two non-resonant production mechanisms: Dodelson- 
Widrow (DW) (sterile-active mixing) and Boson-decay (BD) near the electroweak scale. The 
former yields a distribution function proportional to a thermal fermion but with a decoupling 
temperature Td ~ 150 MeV, whereas the latter leads to a strongly non-thermal distribution 
with a decoupling temperature Td ~ 100 GeV that favors small momentum and yields a 
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colder species of sterile neutrinos for a given mass. For a sterile neutrino with mass ~ keV 
the (DW)-species is warmer with k^ ~ 7.7 (Mpc) -1 and the (BD)-species is colder with 

k^ D ^ ~ 14.12 (Mpc) -1 and its transfer function features a longer plateau for k < kf s as a 
consequence of the ISW enhancement during stage I. 

Although the power spectra is strongly suppressed by free streaming at the scales at which 
(WDM) acoustic oscillations emerge, we conjecture that non-linear gravitational collapse 
may amplify these oscillations into peaks and troughs in the matter distribution at small 
scales, leading to dumpiness on mass scales ~ 1O 9 M for (DW) and ~ 1O 8 M for (BD). 
Perhaps coincident ally this latter scale is of the order of the mass contained within a half-light 
radius in the (DM) halos of spiral, low surface brightness and dwarf spheroidal galaxies [84J. 

Our study also reveals a quasi- degeneracy between the mass, properties of the distribution 
function and decoupling temperature of the (WDM) candidate: particles with the same 
mass but that decoupled at different temperature with very different distribution functions 
may yield similar power spectra in a wide range of scales. As an example of this (quasi) 
degeneracy, the power spectra of (DW) sterile neutrinos with m ~ 2 keV is similar to that of 
a (BD) sterile neutrino with m ~ 1 keV for k < 12 — 15 (Mpc) -1 . This result suggests caveats 
on the constraints on the mass of sterile neutrinos from current (WDM) N-body simulations 
and Lyman-a forest data that typically input the distribution functions of thermal or (DW) 
species. 

We have compared the results for the transfer function for sterile neutrinos produced 
via the (DW) mechanism from the semi-analytic formulation presented here to the results 



obtained in refs.|31l. 1401 . 14 1| from the Boltzmann codes. Although we recognized several 
caveats in the comparison, we find excellent agreement to < 5% between the results from 
the Born approximation ( 1V"I.38j) and the the numerical fit to the result of Boltzmann codes 
presented in ref. 31] in the region of scales where the fit is valid. 

The next step of the program will explore a numerical solution of the full Gilbert equation 
(1VI.3j) along with its comparison to the Born approximation and will be reported elsewhere. 
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Appendix A: Simplification of I\ 

It is convenient to introduce the variables 



x 



krj 4a/2 k 

71 ; ro( " ) = 7r^ >>1 



and change integration variable from u' to r] using eqns. ( IIH.32IIH.8I) ). yielding 



h = -64 

r 1 f Xeq t( \ d f sirl & 

®y Jx NR ^ v x 
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dx 
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where 



U = u + \\n[w{k)\ ; f( x ) = l + ^— 
I w\k) 
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xnr 



Integrating by parts and neglecting terms oc l/w(k) C 1 we find 

l r . sinfx, 
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a v( u ~ 2 ln ( a; )) 



(A.6) 



In the second and third line in the above expression we have approximated f(x) ~ 1 since 
XNn/tt(k) ~ (Kg) 5 <C 1 and the contribution from the upper limit to the integral in the 
third line (the region where x/w{k) ~ 1 ) is suppressed by ~ l/x 2 eq ~ k 2 q /k 2 . It is convenient 
to extract the singular term ocl/xasx~£Ar,R when xnr «C 1 (this is the CDM limit), 
integrating by parts again, leading to 
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(A.7) 



The CDM limit corresponds to a — )■ 0, x NR — )■ 



Appendix B: Numerical implementation of the Born approximation. 



The first step in the numerical implementation is to obtain y y 2 for the given distribution 
function of decoupled WDM particles and to input this value into the mode functions Q, P 
given by eqns. f lV1.261IVL29|) . 

For numerical implementation, it is convenient to take the wavevector k in units of 
(Mpc) -1 and to write 

a = Cl k ; Cl = 0.22f-j ) (B.l) 



k = c 2 k ; c 2 = CiVy (B.2) 
c 22 = ^= = 0.289 c 2 (B.3) 
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along with eqns. f lHI.29IIHI.47p ). lead to 

I r sin(c 22 £;) 
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The last integral term is < 10 3 for k > 0.2 and can be neglected for small scales. 
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